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^ . Abstract 

' ^ I Recent evidence in rodent cerebral cortex and olfactory bulb suggests that short-term dynamics 

of excitatory synaptic transmission is correlated to the occurrence of stereotypical connectivity 
motifs. In particular, it was observed that neurons with short-term facilitating synapses form 
^ ' predominantly reciprocal pairwise connections, while neurons with short-term depressing synapses 

^ ; 

OO ' form unidirectional pairwise connections. The cause of these structural differences in excitatory 

T— I , 

' synaptic microcircuits is unknown. 

I We propose that these connectivity motifs emerge from the interactions between short-term 

synaptic dynamics (SD) and long-term spike-timing dependent plasticity (STDP). While the im- 



X 



pact of STDP on SD was shown in simultaneous neuronal pair recordings in vitro, the mutual 
interactions between STDP and SD in large networks is still the subject of intense research. Our 
approach combines a SD phenomenological model with a STDP model that captures faithfully 
long-term plasticity dependence on both spike times and frequency. As a proof of concept, we 
explore in silico the impact of SD-STDP in recurrent networks of spiking neurons with random ini- 
tial connection efficacies and where synapses are either all short-term facilitating or all depressing. 
For identical background inputs, and as a direct consequence of internally generated activity, we 
find that networks with depressing synapses evolve unidirectional connectivity motifs, while net- 
works with facilitating synapses evolve reciprocal connectivity motifs. The same results hold for 
heterogeneous networks including both facilitating and depressing synapses. Our study highlights 
the conditions under which SD-STDP explains the correlation between facilitation and reciprocal 
connectivity motifs, as well as between depression and unidirectional motifs. These conditions may 
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lead to the design of experiments for the vahdation of the proposed meehanism. 

Author Summary 

Understanding the brain requires unveihng its synaptie wiring diagram, which encodes memo- 
ries and experiences. Synapses are, however, more than mere plugs connecting neurons; they are 
communication channels implemented by biophysical systems that display history-dependent prop- 
erties. During prolonged repeated activation, synapses may functionally undergo "fatigue" or may 
"warm up" . These effects, known as short-term plasticity, alter temporarily the communication 
properties of a synapse, i.e., over tens to hundreds of milliseconds. At the same time, synapses 
undergo structural changes that persist from hours to days, known as long-term plasticity. 

Recent progress in electrophysiology enabled simultaneous access to the synaptic wiring dia- 
gram in microcircuits, and to its short-term features, i.e., the "fatigue" or "warm-up" properties. 
Non-random connectivity patterns were reported, as i) co-occurrence of "fatigue" features and 
unidirectional connections between two neurons, or ii) co-occurrence of "warm-up" features and 
reciprocal connections. We present a biologically plausible explanation for those patterns, based on 
the interaction of short- and long-term synaptic plasticity. We further formulate a computational 
model and provide an interpretation of the simulation results by means of its statistical description. 

Introduction 

Among the most exciting challenges in Neuroscience, the investigation of the brain wiring dia- 
gram known as connectomics made spectacular progresses and generated excitement for its per- 
spectives ([T|). Novel discoveries in molecular biology ([11131111), neuroanatomical methods, ([51 
ini), electrophysiology ([3 [HI O, and imaging ([TUl [TTl [T^ have pushed forward the technological 
limits, for an ultimate access to neuronal connectivity. This is in fact the most important level 
of organization of the brain (|13[) . pivotal to understanding the richness of its high-level cognitive, 
computational, and adaptive properties, as well as its dysfunctions. 

At the microcircuit level ([H [HI [HI [13 [HI), the non-random features of cortical connec- 
tivity have recently raised a lot of interest ([2 The occurrence of stereotypical connectivity 
motifs was experimentally demonstrated and, in some cases, accompanied by physiological infor- 
mation on neuronal and synaptic properties ([3 [IHl HOI HI) : on activity-dependent short-term ([^ 
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[22) and long-term plasticity (|23l ?I^ and rewiring (|25l [25)) . These physiological details are com- 
plementary to anatomical connectivity mapping, and are known to underlie structural, dynamical, 
and computational network properties (I27p . In fact, information transmission between neurons in 
the brain takes place by means of more than mere "connectors" . For instance, short-term dynam- 
ics (SD) of synaptic efficacy is very common and has been quantified as transient and reversible 
facilitation or depression of postsynaptic responses, upon repeated presynaptic activation (j28p . 
Relaying information to downstream neurons in a microcircuit is thus dependent on past-history, 
determined by presynaptic activation frequencies, and shaped by the SD profile at each synapse 

(ED. 

Wc address recent experimental findings obtained in rodent cortices p9p , where short-term facil- 
itation and depression were found to correlate to specific, pairwise, connectivity motifs: neurons, 
connected by synapses exhibiting short-term facilitation, form predominantly reciprocal motifs; 
neurons, connected by synapses exhibiting short-term depression, form unidirectional motifs. In- 
terestingly, the same over expression of connectivity motifs has been observed in another brain 
area, i.e. the excitatory microcircuitry of the olfactory bulb (|29p . The cause of these structural 
differences are largely unknown. 

Inspired by the theory by Clopath et al. (2010) on the relationship between neural code and 
cortical connectivity pop . we hypothesise that interactions between short-term and long-term 
synaptic plasticity could explain the observed pairwise connectivity motifs. In pop . spike-timing 
dependent long-term synaptic plasticity (STDP) was shown to account in silico for the emergence of 
non-random connectivity motifs in networks of model neurons pop . Adopting the same framework, 
we add to it a standard phenomenological description of SD (j5T|) , and study SD-STDP interactions 
in recurrent networks of Integrate-and-Fire neurons p2p . 

Extensive computer simulations of the evolution of synaptic efficacies under both external and 
self-generated activity have revealed an over-expression of reciprocal motifs on facilitating synapses, 
and of unidirectional motifs on depressing synapses. The departure from the initial random con- 
nection efficacies results from the interplay between neuronal activity and synaptic plasticity mech- 
anisms, over long (STDP) and short time scales (SD). Supported by a mean-field analysis ([551 [5^ 
[551 [551 [57]), we conclude that the SD-STDP interplay is controlled by the distribution of firing 
rates in the network, through the switch of STDP from a correlational "pre-post" temporal mode 
at low firing rates, to a "Hebbian" rate mode at high firing rates, earlier reported experimentally 
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by Sjostrom et al. (|38p. Short-term facilitation provides a positive feedback in recurrent niicrocir- 
cuits and it boosts firing, while depression acts as a negative feedback and it effectively attenuates 
high firing rates and discouraging their reverberations. In the "Hebbian" mode, STDP mostly 
gives rise to long-term potentiation (LTP) of synaptic efficacy, which becomes stronger resulting in 
reciprocal motifs. This leads to increasing firing rates of bidirectionally coupled neurons, with LTP 
acting as a long-term positive feedback, until the connection efficacies arc ultimately stabilized by 
saturation. In the "pre-post" (temporal) mode instead, it is the spike timing that drives LTP, 
which strongly favors the emergence of spatially asymmetric, i.e., unidirectional, motifs (I30p . 

Materials and Methods 

We study the properties of networks of excitatory spiking neurons ([5^ . connected via plastic, 
current-based, svnapses ipTl 1511 [551 1551 [501 157)) via a set of simulations. We introduce a convention 
for distinguishing between "strong" and "weak" connections consistent with pop . and a simple 
measure to quantify the occurrence of pair-wise motifs in our simulations. We finally provide a 
Wilson-Cowan firing rate model that is helpful for the interpretation of the numerical results. The 
numerical values of all parameters are indicated in Table [T] 

Neuron model 

The network is composed of identical adaptive exponential Integrate-and-Fire neurons (|32p . each 
described by a membrane potential Vm(t) and by a spike-frequency adaptation variable a;(t) (l39p . 
Below a threshold Ve, Vm{t) satisfies the charge-balance equation 

Cm Vm = gieak (Eleak - Kn) + Qleak Ay 6*^^" ~ ^ - X + Igyn + lext (1) 

where Igyn is the synaptic input from other neurons and lext the external (background) input 
currents. When a spike occurs, i.e., Vm{t) crosses Ve, Vm is reset to a Ereset- 
The spike- frequency adaptation variable x{t) evolves as 

X = a (Vm - Eleak) " X (2) 

When a spike occurs, x evolves as a; — ?• x + A^. The numerical integration of Eq. |T]is suspended 
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for a period of time Tarp following each spike, to mimic absolute refractoriness, during which Vm 
remains "clamped" at Ereset- 

The model details are not essential to our conclusions (see Supplemental Information, Figure 
S7). 

External (background) inputs 

Each neuron is identified by an index i = 1,2,3... and it receives a time variant input lext i 
according to the following protocols. 

Toy Network: ( e.g., Fig.[Tj3, 10 neurons) lext i consists of a 0.5 nA constant current, as well 
as periodic 1 nA square pulses. Pulses occur as in a traveling wave of activity, which moves every 
5 msec from one unit, e.g., the i-th neuron, to its next index neighbour, e.g., the (i + l)-th neuron, 
see also (j30p . Each pulse is delivered in turn to all neuronal indexes, as an extremely narrow bell- 
shaped profile with unitary amplitude and standard-deviation of 0.5, resulting in neighbouring 
neurons being only slightly stimulated. Each pulse is of sufficient amplitude to elicit neuronal 
firing in the unit where the bell-shaped profile is centred on, e.g., the i-th unit (unless refractory). 

Large Network: (e.g.. Fig. [SJ 1000 neurons) I^xt is as in the toy network with the addition of 
an uncorrelated gaussian noisy current (|40p. with mean /i, standard deviation a = 200 pA, and 
autocorrelation time length r^yn ~ 5 msec. Parameter fi is drawn randomly before launching the 
simulation and for each neuron of the network, from a gaussian distribution with mean 200 pA 
and unitary coefficient of variation. The noisy current mimics asynchronous synaptic inputs from 
(not explicitly modelled) background populations (|in W^ . 

Internal (synaptic) inputs 

Neurons connect to each other according to a fixed anatomical wiring matrix [Cy], which indicates 
whether neuron j-th projects to neuron i-th, i.e., C'ij = 1, or not, i.e., Cy = 0. The matrix [C'ij] is 
obtained from an all-to-all connectivity without autapses (i.e., Cu = 0), upon randomly pruning 
~ 20% of its elements (see e.g. Fig. [TJ3), to introduce variability in neuronal firing. A more 
substantial reduction of the structural connectivity does not affect qualitatively our conclusion, 
although it downscales the number of plastic synapses available for statistical analysis. This would 
require larger networks to be simulated in order to obtain the same statistical confidence over our 
results, and would require upscaling Gij (e.g., through increasing Wmax, see the next subsections) 
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to obtain the same network firing rates. 

The z-th neuron receives at any time a synaptic current I i, described as 

N CO 

^syn i — ^syn i / '^syn ~1~ ^ ^ ^ ^ '^(^ ) (^) 

i=i / 

where represents the occurrence time of the f-th spike emitted by the j-th presynaptic neuron, 
and Gij the amphtude of the postsynaptic current (PSC), corresponding to the activation of the 
synapse by the presynaptic j-th neuron. The Dirac's delta function S{t) represents the occurrence 
of a presynaptic action potential. Eq. [3] models incoming individual PSCs as waveforms with 
instantaneous rise time and slower decay time (j43|) . imposes a linear postsynaptic superposition of 
effects, and implies that in the lack of any presynaptic spiking activity, Isyn i decays exponentially 
to zero with a time-constant Tsyn, and when a presynaptic spike is fired it evolves as Isyn i 

Frequency-dependent short-term synaptic dynamics (SD) 

Gij defines the amplitude of the PSC from presynaptic neuron j-th to postsynaptic neuron i-th. 
On short timescales, its value changes as a function of the activation history of the presynaptic 
neuron, leading to transient and reversible depression or facilitation of synaptic efficacy (j2ip . This 
is referred here as homosynaptic plasticity and implemented by having Gij proportional to the 
amount of used resources for neurotransmission Uij Vij and their maximal availability Aij, i.e., 

Gij -^ij '^'ij ^ij • 

Frequency-dependent short-term synaptic dynamics is described by the following equations, 
with a different set of parameter values to capture depressing or facilitating synapses (pijl: 



i^-r^j)/Trec - Y,'^^ U^-j Tij 5{t ~ tk,) 

(4) 

-U^-j/Tfacil + U (1 ~ Uij) S{t -tk^) 



For the sake of notation, indexes have been dropped from parameter U as well as from the time 
constants Tree and Tfaeii in Eqs. 21 although in general each synapse has its own parameters (see 
Tabled]). Eqs. S] reduce to the following update rules: i) when no spike is fired by the presynaptic 
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neuron j, Uij and recover exponentially to their resting values, U and 1. respectively; ii) 
as a presynaptic spike occurs, r^j is reduced as — s> (1 — Uij) rij, while u^j is increased as 
Uij — Uij + Uij (1 — Uij). The impact of short-term plasticity of PSCs amplitude is exemplified 
in Fig.fflCF. 

Spike-timing dependent long-term plasticity (STDP) 

We further extend the description of PSCs (Eqs. by an additional scaling factor Wij. to 

incorporate STDP ([23]), see also 

Gij = Wij Aij Uij Tij. 

(5) 

Wij changes on timescales longer than r^ec and Tfadi and according to the correlated activity of 
both pre- and postsynaptic neurons, as in capturing spike-triplets effects ([551 [5(11 [57)) . This 
is referred here as heterosynaptic plasticity. Each neuron is complemented by four variables, i.e., 
91 ; 92 ; oi , 02, which act as running estimates of the neuron firing rate, averaged over distinct time 
scales, i.e., t^^, t^^, Tqi, To2- In the lack of any firing activity of the j-th neuron, those variables 
exponentially relax to zero: 

while each time the neuron fires, these variables are increased by a unit: 

91, 91, + 1 92, ^ 92, +1 oi^ ^ oi^. +1 02, 02, +1. (7) 

This enables a compact formulation of STDP: when the j-th neuron spikes, the following updates 
are performed for all the indexes i: 

Wy ^ W^j -T] oi^{t) [A^ + Ajq2^{t~e)] 

(8) 

Wj^ ^ Wj, +T] qi^(t) [A+ + A+02,(t-e)] 

as the j-th neuron is presynaptic to all connected i-ih neurons, and postsynaptic to all connected 
j-th neurons, respectively. Numerically, the evaluation of 92, and 02, is performed just before the 
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neuron j-th spikes, as indicated by the infinitesimal time-advance notation of e. When no spike 
occurs, Wij maintains indefinitely its value. The instantaneous value of each Wij is bounded in the 
range [0 ; Wmax]- Unless otherwise stated, Wij is initialized from a uniform random distribution 
between and Wmax, prior to the start of each simulation. 

We remark that our choice for combining STDP and SD together aims at a phenomenological 
description of their interactions. It is assumed that all STDP variables are local as well as that 
STDP and SD models act independently. This allows us to easily analyse network interactions and 
avoid the paradox of pre-synaptic firing detection under strong short-term depression of synaptic 
efficacy. In fact, we assume for simplicity that the STDP coincidence-detector is located at each 
synapse and that such a detector is able to access timing information of presynaptic spikes, re- 
gardless on whether they result in large or small PSCs. Similarly, the same coincidence-detector 
will also sense all the spikes emitted by the post-synaptic neuron. This is the only information 
required to calculate the long-lasting synaptic modification. 

Convention on "strong" and "weak" connections and network symmetry index 

In this study, we focus on the appearance or disappearance of a strong connection between two 
neurons, only for units that are anatomically connected, i.e., Cy = 1. For the sake of comparison, 
we adopted the framework of Clopath et al. (2010), where the activity-dependent appearance or 
disappearance of a connection conventionally occurs in terms of a competition among the "strong" 
links in a "sea" of weak synapses ([7]). As in their paper, we adopt the convention of identifying as 
"strong" those connections whose corresponding STDP synaptic factor Wij is above the 2/3 of its 
upper bound W^ax- 

With such a definition, we measure the symmetry in the network connectivity by counting 
reciprocal or unidirectional motifs as (but see e.g., (j45|) for alternative definitions) 



where N is the size of the matrix W , as well as the size of the network. The symmetry index siW) 
takes values in the range [0 ; 1] and depends on the average similarity between elements of W that 
are on symmetric positions with respect to the diagonal. Following our previous convention, Wij 
and Wji are first normalized and then zero-clipped: Wij = Wij / Wmax if Wij > 2/3 Wmax, 




(9) 



i=l 
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and otherwise W*j ~ 0. In Eg. IS46[ M represents the number of null pairs {W*p Wj^} — 
{0, 0} that occur as a consequence of clipping or by initialization. Evaluating s{W) on networks 
with a majority of unidirectional connections results in values close to (e.g., Fig. while 
its evaluation on networks with a majority of reciprocal connections results in values close to 1 
(e.g., Fig. [2j3). For uniform random matrices W, the full statistics of s{W) can be calculated 
(see Supplemental Information) and employed for deriving a significance measure for s as a p- 
value, representing the probability that the value of s observed in simulations could result by 
chance. This symmetry measure is useful only in silico, where the entire connectivity matrix 
as well as its maximal synaptic efficacy are known. The same measure cannot be considered 
as an operational tool for quantifying experimentally-extracted microcircuit connectivity, where 
statistical counting of connectivity motifs and their comparison to chance- level are relevant ([7]). 
Approaching a qualitative comparison to the experiments of Wang et al. (2006), such a statistical 
counting has been used in the large simulations of Fig. [5l 

Mean-field Network Description 

We approximate the firing rate of the network of Integratc-and-Fire units through its mean-field 
dynamical description (|^ Wf\ IM)) . closely following earlier work (PTl We consider the sim- 
plifying hypotheses that i) the network consists of one or more non overlapping subpopulations of 
excitatory neurons (see Figs. |3K-B and Fig. [SJfV) and that ii) neurons within each subpopulation 
share identical synaptic coupling, connectivity, and short-term synaptic plasticity properties, i.e., 
all depressing or all facilitating, as in Fig.[51A_; see Fig.[71A. for an exception. We also assume that iii) 
for each of the (sub)populations, the individual neuronal firing follows a Poisson distribution, with 
instantaneous mean firing rate E(h), which depends monotonically on the corresponding overall 
input currents h. Under these hypotheses, neurons can be only distinguished by the subpopulation 
they belong to, i.e., depressing D or facilitating F, and their firing can be identified as ED{hi)) 
and Ep{hF). For the case of two subpopulations (Fig.lHK), Hd and hp evolve over a characteristic 
time scale t as 



T ho 



ho + Jdd En + Jdf Ep + I^xt 



< 



(10) 



T hp 



hp + Jpp, Ep) + Jpp Ep + lext 
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where lext represents the external input, Jdd and Jpp the average synaptic efficacies of recurrent 
connections within each subpopulation, and J^p ^^nd Jpp) the average synaptic efficacies of inter- 
subpopulations connections. On a first approximation, J can be considered as the ensemble average 
of Gij. Firing rates Ed and Ep are computed from ho and hp as threshold- linear frequency- 
current response functions: Ep, = [hpi — 9]^ and Ep = [hp — 0]^, with [x]^ = max{x;0} 
(for alternatives see (|48l |49)) ). Jdd, Jff, Jdf, Jfd undergo plastic changes, on both the short 
and long time scales: indicating by Jab the coupling between the presynaptic population b and 
the postsynaptic population a, then Jab ~ Wab Aab Ub Xb with a,b £ {D,F}. The mean-field 
variables Ub and Xb depend only on the presynaptic firing rate Eb, and capture the short-term 
homosynaptic plasticity in the mean-field approximation of Egs. l4l ( j3ip . as: 

UD = {Ud - Ud) /Tfacti D + Ud (l - ud) Ed 

Xd = {I - Xd) /Tree D + UD Xd Ed 

(11) 

Up = {Up - Up) /Tfaeii F + Up {I - Ud) Ed 

Xp = {1 — Xp) /Tree F + Up Xp Ep 

On longer timescales, the mean-field approximation of STDP given in (I35p is adopted here by 
altering the factor Wab, which evolves as a function of both presynaptic Eb and postsynaptic Ea 
firing rates: 

i Wab = -A^ To, Eb Ea - To, Tg, E^ Ea + A+ T,, Eb Ea + A+ T,, T^, Eb E^ . (12) 

Results 

We hypothesize that synaptic connectivity refiects the interaction of short-term plasticity, such 
as homosynaptic depression and facilitation (j2ip . with long-term plasticity. Inspired by a recent 
theory, which links connectivity to spiking activity and to the neural code (j30p , we ask the question: 
could known relationships between short-term plasticity and spiking activity in recurrent networks 
(|5n [551 explain the occurrence of connectivity motifs (IT^ observed experimentally? In 



Connectivity Motifs by Synaptic Plasticity 



11 



support to our hypothesis, we present simulations of networks of excitatory neuron models with 
activity- dependent plastic connections. We then study the necessary conditions for the emergence 
of these specific connectivity motifs by standard firing rate models. 

A toy microcircuit model 

In order to demonstrate and test the key features of our hypothesis, we consider a simplified rep- 
resentation of excitatory neuronal microcircuits, as a network of adaptive exponential Integrate- 
and-Fire units ((5^ . The validity of our results does not depend on the specific details of the model 
neuron and of its spike-frequency adaptation parameters chosen here (Supplemental Information, 
Figure S7). Neurons are connected to each other through excitatory synapses (Fig. [T)3), whose 
efficacy undergoes either short- or long-term plasticity. Using a standard model for short-term 
synaptic dynamics (SD) (j2ip . we simulate homosynaptic, use-dependent modifications of the post- 
synaptic potentials (PSPs) amplitude. We also employ a recently introduced phenomenological 
description of heterosynaptic long-lasting potentiation or depression of PSPs, able to capture with 
great accuracy both spike-timing and frequency effects ([551) . We refer to this long-term associative 
plasticity as spike-timing dependent (STDP). Both SD and STDP occur simultaneously with neu- 
ronal dynamics, although across distinct timescales. We note that STDP interacts with short-term 
plasticity as a frequency-independent scaling factor of PSPs amplitude, rather than contributing 
to a redistribution of synaptic efficacy in the sense of (jSOl [5T|) (see Discussion). In addition to 
excitatory PSPs evoked by trains of presynaptic spikes, neurons receive an external current input, 
deterministically played back over and over, as a traveling wave of activity (Fig. [1)3). Such an exter- 
nal stimulation recreates the temporal coding of Clopath et al. (2010), which imposes deterministic 
spike-timing correlations among evoked spikes (see Methods) . The external input can be regarded 
an oversimplified generic e.g., thalamic, input with known temporally correlated structure. 

We define two microcircuits, identical for all aspects of neuronal properties, maximal synaptic 
efficacy, anatomical connectivity, and external inputs, with the exception of the SD properties. 
Specifically, one microcircuit includes exclusively short-term facilitating synapses (Fig. [Tp), while 
the other includes only depressing synapses (Fig. [Tf). In order to describe changes in the micro- 
circuit connectivity, we adopt the same framework of Clopath et al. (2010) where connections 
form or disappear via competition among the "strong" links in a "sea" of weak synapses ([7|) (see 
Methods) . The connectivity, which is randomly initialized, slowly evolves into largely non-random 
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configurations during the simulation (Fig. [Tp,G). At the steady-state, these configurations match 
the experimentally observed correlations: reciprocal motifs emerge in cell pairs more often than 
unidirectional motifs, when synapses are facilitating; the opposite occurs when synapses are de- 
pressing. This is revealed both by direct inspection of the synaptic efhcacy matrix [Wij] (e.g., see 
Fig. |8]) and by quantification of its symmetry index s (see Methods). Although other statistical 
measures can be defined, as for instance the percentage of weak elements in [Wij] or the sparseness 
of the emerging connectivity, here we focus on the reciprocal versus unidirectional features of its 
strong elements, according to the same convention of Clopath et al. (2010). When s takes values 
close to 1, almost all of the existing pairwise anatomical connections are reciprocal and \Wij\ is 
almost a symmetric matrix. On the other hand, for values of s close to 0, unidirectional or very 
weak connection motifs prevail. 

Wc underline that the prevalence of external inputs over recurrent synaptic inputs (or vice 
versa) is not imposed a priori but it depends on the interaction between SD and STDP. For 
facilitating microcircuits, the overall recurrent (synaptic) input that a neuron receives becomes 
progressively larger than the external input, as soon as more and more reciprocal connectivity 
motifs are established by STDP. For depressing microcircuits, the opposite holds, as more and 
more connections are weakened by STDP. 

These results, summarized in Fig. [T] for a sample microcircuit composed by ten neurons, 
have been confirmed over 2000 repeated simulations, analysed in Fig. [51 Over a slow timescale, 
STDP results in a very small degree of symmetry (s = 0.01 ± 0.01, mean ± stdcv), with 
high significance (jp < 10^^) for each of the repeated simulations involving depressing synap- 
tic connections (Fig. HP). Under identical external inputs to the microcircuit, STDP leads in- 
stead to a large degree of symmetry (s — 0.61 ± 0.10, mean ± stdev), with high signif- 
icance {p < 10"**) in about 75% of the simulations involving facilitating synaptic connec- 
tions (Fig. HP). In the remaining 25% of the cases, the resulting symmetry value was in the 
range [0.25; 0.55], suggesting that to some extent, the variability observed in experiments (dHl 
HHl) was replicated in the simulations, where not 100% of the times the motifs correlations emerged. 
In this simplified example, the variability is attributed mostly to the sparse anatomical connectiv- 
ity and the irregular firing regimes this imposes. Doubling the simulation time (not shown) led to 
a very minor reduction of the variability on s, by less than 3% and only for the facilitating synaptic 
connections, suggesting that stationarity had been already reached. 
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A simple mechanism for the emergence of motifs 

In the microcircuits considered above, the statistics of the anatomical connectivity, the external 
inputs, as well as the single-neuron properties are identical. Any difference in long-term modifica- 
tion of connections could only arise i) from differences in the collective neuronal activity and take 
place ii) through the spike-timing and firing-rate dependent mechanisms underlying connection 
motifs formation, as first described in pop . In the microcircuits considered here, neuronal activ- 
ity is known to depend on the connectivity and the short-term changes in PSPs amplitude pip . 
Indeed, during the 2000 simulations of Fig. [2J the microcircuits including short-term depressing 
synapses collectively fired at 20 ± Hz (see also Fig. [Tfl), lower than the recurrent microcir- 
cuits employing short-term facilitating synapses, which fired at 59 ± 4 Hz (see also Fig. HJi;). 
Facilitating synapses charge up with ongoing presynaptic activity and they recruit more connected 
neurons, as in a positive feedback. On the contrary, depressing synapses get soon fatigued, re- 
sulting into a weaker subsequent recruitment of postsynaptic neurons, whose firing would further 
depress synaptic efficacy, as in a negative feedback. 

We identify such an asymmetry as the cause underlying motifs emergence, and the firing-rate 
dependence of STDP (|24)) as its mechanism. In particular, it is the functional switch of STDP from 
a correlational "pre-post" temporal mode at low firing rates, to a "Hebbian" rate mode at high 
firing rates, that entirely accounts for the asymmetry in the connectivity motifs. In addition to the 
spike-timing information (Fig.[2K), the STDP model employed here faithfully captures the strong 
dependency on the neuronal firing rates (Fig. 03). For the sake of illustration, we isolate the 
impact on synaptic efficacy of Eas. l6l7l51 in a two-neuron system, with one neuron projecting to the 
other via a single synapse. We simulate the long-term change in PSPs amplitude at that synapse 
(sec Fig.[5}3 and psp ). imposing 75 pre-post spike pairing events, evoked at different uniform firing 
frequencies. We study two cases: (i) each presynaptic spike precedes the postsynaptic spike by 
10 msec, i.e., tpre < tpost, or (ii) vice versa, i.e., tpre > tpost- In agreement with the experiments 
(PS)) . above 30 — 40 Hz long-term potentiation (LTP) prevails on long-term depression (LTD), even 
in those cases when spike-timing per se would promote LTD. i.e., tpre > tpost- Below that critical 
frequency, LTP or LTD refiects causal or anti-causal relationships between pre- and post-synaptic 
firing times, respectively p4p . 

It follows that facilitating microcircuits would preferentially develop reciprocal connectivity 
motifs, irrespectively of the spike-timing, due to their higher firing rates that promote LTP. The 
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consequent increased recurrence in the connections would further increase the overaU population 
firing rate, above the critical frequency. This increment in the firing rates leads at first place to 
LTP. Depressing microcircuits instead do not develop preferentially reciprocal motifs, but reflect 
the asymmetric temporal structure of the input. On a first approximation, their connections 
would rather not promote persistent high activity regimes; instead they mostly respond to external 
stimuli. These stimuli are imposed here at low rates and give rise to feed-forward links, due to the 
repeating traveling- wave features of the stimulation, as in pop . 

These specific connectivity and activity configurations arc stable when imposing reflecting 
boundaries, as in the numerical implementation of STDP ([521 [STl |53)) . These boundaries limit 
the synaptic efficacy to a minimal and a maximal values (see Methods). We also remark that the 
firing rate dependence of STDP, shown in Fig. [2j3, is not captured by all the models proposed in 
the literature. In our case, the minimal set of STDP model features sufficient for motifs emergence 
must include the reversal of LTD into LTP at high firing rates. 

As a conclusive illustration of the simple mechanism for motifs emergence, we performed addi- 
tional negative and positive control simulations (Fig. [3]) , studying the impact of alternative formu- 
lations of STDP. By setting A'^ = = 0, and A^ = 4.5 10"^ in Eq.[5]and slightly modifying 
Eq. [7] (see the Supplemental Information) one can reproduce the pair-based STDP model ((5^ 
I54p . When parameters are adjusted so that this model has an identical spike-timing dependency 
(i.e., compare Fig.[3J\ to Fig.l^K), the resulting frequency dependency is completely different (i.e., 
compare Fig. to Fig. 03), showing no reversal of LTD into LTP at high frequencies. Vice 
versa, by inverting signs and swapping the values A2 , A^ , A2 , and A^ in the same equations, 
it is possible to ad-hoc reverse the temporal dependency of STDP, as described experimentally 
for anti-STDP (aSTDP) ([55]) . while leaving the frequency dependence roughly intact (i.e., com- 
pare Fig. 1313, D to Fig. [5]\,B) and featuring the reversal of LTD into LTP at high frequencies. Of 
course, this modified "triplet rule" model should not be generalized as novel accurate description 
of aSTDP. since more data would be anyway needed to access its firing rate dependence. Panels 
E-H of Fig. [3] repeat the simulations of Fig. [5J and demonstrate that only when the frequency- 
dependence of STDP is realistic ([35]), the heterogeneity in the network firing rates leads to the 
emergence of asymmetric connectivity motifs (compare Fig. [3fl,G or Fig.[3p,H to Fig. [2P,D). 
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Large homogeneous microcircuits 

Due to the simplicity of the mechanism disclosed above, we can analyse larger populations of 
neurons interacting via homogeneous and heterogeneous short-term plastic synaptic connections. 
Prior to running numerical simulations, we qualitatively investigate the generality of our results 
and the conditions for the asymmetry in the emerging firing rates. We carry out this analysis 
by a standard Wilson-Cowan firing rate description of neuronal networks dynamics (see Eq. llOp 
(|Tf)) . This approach has been already extended to the case of short-term synaptic depression 
and facilitation (sec Eq. [TTI) ([331 [Ml IMl), long-term Hcbbian plasticity (gT]), or both (|SH|). This 
technique can be used as long as the characteristic times of long-lasting plasticity are much longer 
than those of neural and synaptic short-term dynamics, as in our case. 

For the sake of simplicity, we initially omit STDP (Eq. and consider recurrent networks 
of neurons connected by homogeneous synapses, facilitating (Fig.|4]A.) or depressing (Fig.|4j3). We 
sketch these networks as oriented graphs, where arrows represent connections with average efficacy 
JpF or Jdd, correspondingly. We do not explicitly include inhibition, which, nevertheless, does 
not qualitatively alter the validity of our results, as examined in the Supplemental Information. 
Instead, we distinguish whether the average external input to a generic neuron is zero, i.e., referred 
to as balanced inputs (lext ~ in Eq. llOp . or is set to a positive value, i.e., referred to as 
unbalanced inputs, (lext =^ 5 in Eq. llOp. The population firing rates can be studied via standard 
methods as in dynamical system theory (|59p . i.e., analyzing the system of Eq. 1101 and evaluating 
their steady-state solutions. 

For instance, given a specific combination of external inputs and recurrent average synaptic 
efficacies, i.e., Jff = Jdd = 4 for unbalanced inputs and Jff = Jdd ~ 10 for balanced 
inputs, the steady-state solutions of Eq. [TU] provides the equilibrium firing rates. These can be 
then compared to the critical frequency of STDP (see Fig. [2}3), invoking the same mechanisms for 
motifs emergence proposed by Clopath et ai, (2010): if above the critical frequency, then LTP 
prevails and STDP will reinforce reciprocal synaptic efficacy; if below the critical frequency, then 
LTP or LTD will be determined by spike-timing information. Panels C and D in Fig. 2] illustrate 
graphically in the plane E, h, the actual location of the steady-state solutions of Eq. (TU] These 
are derived by means of geometrical analysis techniques, similarly to the nuUclines intersection 
methods (|59p . In fact, at the equilibrium Eq. [TU]can be rearranged so that its roots correspond 
to the intersections between the unitary slope line (dashed) and a given curve (continuous lines; 
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see Supplemental Information). The shape and bending of such a curve depend on the external 
input /e^jt and on the SD parameters, so that the firing rates emerging in the recurrent networks 
of Fig. |3J\-B can be compared to each other (1551 [551 155)) . The dynamical stability of each steady- 
state solution is displayed by a different marker symbols: circles for stable, squares for unstable 
equilibrium points. Similarly to Fig. [2j3, the approximate location of the STDP critical frequency 
has been indicated by a gray shading. Although we did not chose the parameters of the rate models 
to quantitatively match the Integrate-and-Fire simulations (|40l[57)) . we conclude that homogeneous 
facilitating networks generally fire at higher firing rates than depressing networks, for the same set 
of average synaptic coupling and external inputs, and when engaged in reverberating activity (|40p . 

Mathematically the two networks share the same mean-field description, Eqs. [TOlfTTl although 
their numerical parameters are considerably different (see Table [1]) ^T9[ (29)) . Specifically, short- 
term depressing synapses have a much larger recovery time constant from depression (r^ec) than 
facilitating synapses. Facilitating synapses have instead a much larger recovery time constant from 
facilitation {rfacu) than depressing synapses. 

Along these lines and by studying the asymptotic limits of the mean-field equations, we found 
heuristically that the dominating parameter is r^J.- This sets an upper limit to the location of any 
possible firing rate of each network. Then, a network with low values of t^J,, i.e., where the time 
scale of recovery from depression is very long, would fire slower than a network with comparatively 
higher values of T^rJ,, i.e., where the time scale of recovery from depression is very fast or negligible. 
These two cases correspond to the depressing and facilitating networks, respectively (Figs. |3]A_-B). 

These general considerations have been validated and confirmed in networks of 1000 neurons, 
see Methods. In these simulations, neurons are connected with maximal synaptic efficacy A = 6 pA. 
In order to increase the biological realism, in these large scale simulations we introduce fluctuating 
random inputs to each neuron, mimicking background network activity (|42p . Therefore, each neu- 
ron receives an uncorrelated noisy current, as well as a periodic wave-like stimuli (see Methods). 
Figs. (SIA.-B display the count of the occurrence of unidirectional versus reciprocal connectivity 
motifs. Similar to the small toy model, analysed in Fig. [21 unidirectional depressing connec- 
tions significantly outnumber the reciprocal depressing connections, while facilitating reciprocal 
connections prevail on unidirectional facilitating connections. As indicated by Figs. [Sp-E, the 
distributions of the firing rates of the two networks feature the same heterogeneous distributions 
of firing rates: networks of homogeneous depressing short-term plastic synapses fire at generally 
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low rates, while networks of homogeneous faeilitating synapses fire at higher rates. Finally, the 
symmetry index s, eomputed after a very long simulation run, results in a value of 0.28 for the 
depressing network and of 0.99 for the facilitatory network. 

Heterogeneous microcircuits 

We further study the more general case of a heterogeneous network (Fig. [GK) with two subpopu- 
lations: one with faeilitating synapses, and the other with depressing synapses. Arrows represent 
synaptic connections, with average efficacies JfFi Jdd, Jfd, and Jdf- Synapses established 
within neuronal pairs that belong to the same subpopulation, share, by definition, the same SD 
properties, i.e., short-term depressing or short-term facilitating, but not both simultaneously. On 
the contrary, synapses established within neuronal pairs that belong to distinct subpopulations 
have, by definition, heterogeneous SD properties. Thus a total of five categories of connectivity 
motifs are possible in this network: facilitatory reciprocal motifs, depressing reciprocal motifs, 
facilitatory unidirectional motifs, depressing unidirectional motifs, and reciprocal motifs with both 
facilitation and depression. For the first four categories, experiments support strong non-random 
occurrences (|191 . For the case of reciprocal motifs with both facilitation and depression, no 
extensive experimental information has been published. Our results suggest that non-random oc- 
currences of the first four categories arises from SD-STDP interactions, and predict that the last 
category should be largely underexpressed, compared to chance level. 

We first examine the impact of STDP in a simplified two-neuron system, with one neuron 
projecting to the other via a single synapse. Figs. in}3,D, show the long-term change in PSPs 
amplitude as a function of the pre- and postsynaptic firing rates, at that single synapse. Since in 
a heterogeneous network pre- and postsynaptic firing frequencies may differ, we swept the firing 
frequencies of the two neurons throughout all the possible combinations, within a realistic range. 
We study two cases: each presynaptic spike precedes the postsynaptic spike by 10 msec, i.e., tpre < 
tpost, or vice versa, i.e., tpre > tpost- We emphasise that only synapses established within neuronal 
pairs that belong to distinct subpopulations can experience heterogeneous pre- and postsynaptic 
firing rates. In this case however, the impact of spike timing information becomes negligible as 
soon as pre- and postsynaptic neurons fire at different frequencies. In the small minority of cases 
where this is not true, pre- and postsynaptic frequencies are integer (sub)multiples of each other, 
and a transient synchronization of spike times occurs periodically. In these circumstances, the 



Connectivity Motifs by Synaptic Plasticity 



18 



timing information has a specific impact, as revealed graphically by bright or dark lines in the 
plots of Figs.[nj3,D. In all other cases, overall plasticity profiles reflect the conventional associative 
Hebbian LTP/LTD and its consequences ([SOI [CTl [57)1 . 

Intuitively, the heterogeneous population of Fig. |6K can lead to the emergence of connectivity 
motifs. To illustrate this point, we first ignore that subpopulations might interfere with each other's 
firing rate. We assume that the facilitating and depressing subnetworks would be still characterized 
by higher or lower firing rates, respectively, as previously presented for homogeneous networks. 
Then, the LTP/LTD maps of Figs. [6j3,D suggest that such an initial asymmetry of emerging 
firing rates can be maintained indefinitely. The connections Jjjp are increasingly weakened and 
the connections JpD strengthen. The resulting configuration, sketched in Fig. |6lC, is stable. We 
tested and confirmed this statement under the mean-field hypothesis, by studying the dynamics of 
Eqs. 1101 lll[ and 1121 in addition to computing their equilibria. Fig. |B|il-F display the mean firing 
rates of each subnetwork, and the initial time course of synaptic efficacies, with Jjjp progressively 
becoming weaker. This is true only when inter-population efficacies, i.e., Jfd and Jdf, are 
initialized to slightly weaker values than the intra-population efficacies, i.e., Jff and Jdd- Such 
an initial difference between inter-population and intra-population couplings could however emerge 
from fully homogeneous couplings, i.e., Jff = Jdd = Jdf = Jfd = 1, as demonstrated in 
the Supplemental Information. 

We further confirmed that the synaptic configuration, indicated by the mean field models, is 
present in large microscopic numerical simulations, as in Fig.[5lC (see Methods). These simulations 
involve 1000 identical Integrate-and-Fire units, subdivided in two subpopulations of equal size, with 
their anatomical connectivity set to 80% of all possible connections, similar to the homogenous 
case. The maximum synaptic efficacy is set to yl = 12 pA and the initial synaptic connections Wij 
are randomly initialized. As in the mean-field model, the inter-population terms Wij are initialized 
to weaker values than the intra-population terms (but see Supplemental Information for alterna- 
tives) . Each neuron receives an uncorrelatcd background noisy current as well as periodic wave- like 
stimuli, similar to the homogeneous case. As indicated by Fig. [Sf, the distributions of the spiking 
frequencies of individual neurons of the two subnetworks feature the same heterogeneous distribu- 
tion of firing rates as in the previously studied homogeneous networks: subnetworks of depressing 
short-term plastic synapses fire at generally low rates, while subnetworks of facilitating synapses 
fire at higher rates. The location of the critical firing frequency for the STDP is represented again 
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as a grey shaded area. 

Results in Fig. [Sp show all the possible synaptic combinations. As in the data of Wang et 
al. (2006); reciprocal motifs are significantly co-expressed with facilitatory synapses and unidirec- 
tional motifs with depressing synapses. The actual motif count is compared to the null-hypothesis 
of having statistical independence between the connection occurrence within a pair of neurons, 
estimated at a 95% confidence interval upon the same hypothesis of Bernoulli repeated, inde- 
pendent, elementary events. The frequency Q of observing a connection between two neurons, 
regardless of its SD properties, is first estimated by direct inspection of the connectivity matrix 
[Wij]. Then the conditional occurrence frequencies of a facilitatory synapse Qf and of a depress- 
ing synapse Qjj = [1 — Qp) are computed, given that a connection exist between two neurons. 
The null hypothesis for each possible combination is given by standard probability calculus, under 
the hypothesis of independence of the identical events. For instance, the occurrence frequency 
of observing by chance no connections within a neuronal pair is (1 — Q)^, while the occurrence 
of observing by chance a reciprocal motifs with mixed depressing and facilitating properties is 
2 (Q^ Qp Qd)- Finally, in this example, the symmetry index s, computed after a very long simu- 
lation run, resulted in a value of 0.18 for the depressing subnetwork and of 0.66 for the facilitation 
subnetwork. 

Microcircuits with overlapping SD properties 

In the heterogeneous network of Fig. [51 as well as in the homogeneous networks, we make the 
assumption that the SD profile is determined primarily by the identity of the projecting neuron. 
This has been experimentally found in the olfactory, visual, and somatosensory cortices as well as 
in other brain areas (|62l [551 [Ml I55|) . Nonetheless, the assumption on the projection-cell specificity 
can be removed in order to theoretically explore the impact of SD heterogeneity across distinct 
synaptic connections established by the same presynaptic neuron (|66p . 

We assumed that a generic neuron had a certain probability pu ol establishing a short-term 
depressing synapse with a target neuron, and a probability 1 — p^ of establishing a short-term 
facilitating synapse with another one. In this case, individual neurons are still indistinguishable 
and their mean-field synaptic input can be described mathematically Fig. [TlfV (see Supplemental 
Information). 

For small values of pu, the emerging firing rates approximate those of a network of facilitating 
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synapses, while for large values ofpo the firing rates behave as for a network of depressing synapses. 
In other words, the mixed networks behave dynamically as an intermediate case between two 
extremes. This result is quantified in figure Fig. [7J3, where the location of the stable equilibrium 
points has been analysed under the mean- field hypotheses and plotted as a function of pD, for 
different external inputs regimes. The qualitative location of the critical firing frequency for the 
STDP is represented as a grey shaded area. In this situation, and as an explicit consequence of 
the lack of any structure, i.e., compare Fig.[7lA. with Fig.|6]4_, STDP fails to discriminate individual 
connections within the network, but rather shapes them as reciprocal or as unidirectional motifs, 
depending on the particular choice oi po- 
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Discussion 

Our results indicate that time- and frequency-dependent STDP mechanisms may be responsible, 
through internally generated spiking activity in recurrent network architectures, for the observation 
that excitatory neurons connected by short-term facilitating synapses are more likely to form 
reciprocal connections, while neurons connected by short-term depressing synapses are more likely 
to form unidirectional connections. More specifically: 

1. the internally generated firing rates in model networks with facilitating connections arc higher 
than in networks with depressing connections, under identical background inputs; 

2. neurons, participating to such an internally generated activity, are likely to form bidirectional 
connections with each others when firing at sufficiently high rates, reflecting the "Hcbbian" 
mode of STDP; neurons firing at low rates are likely to form unidirectional connections, 
reflecting the temporally asymmetric "pre-post" temporal mode of STDP; 

3. once formed, these connectivity motifs persist and self-sustain themselves through the internal 
firing regimes of the network, from which the motifs emerged from; 

4. externally generated inputs, strongly depolarizing or strongly hyperpolarizing individual neu- 
rons, when prevailing over internally generated activity, may lead to opposite motifs emer- 
gence consistent with Clopath et al. (2010). 

Relationships to previous work and additional mechanisms 

The impact of long-term associative synaptic plasticity in recurrent networks of spiking neurons 
has been earlier studied by many investigators, who proposed mechanisms for the unsupervised 
formation of stimulus-driven dynamical attractors of network activity, in the context of working- 
memory states (j57p . Within the same aims, the interactions between long-term plasticity and SD 
were also already partly explored, both in numerical simulations and in mean-field descriptions 
(155)) . Here, we focused on specific long-term plasticity mechanisms (STDP) (12^ to study the 
emergence of network structure (p71 [551 IS^ 170]) and connectivity motifs ([501 [7T|) . To the best of 
our knowledge, this is the first time that SD is viewed as key element for the emergence of network 
connectivity, and shown to capture, to a certain extent, the experimental observations. 

Our study builds upon assumptions of Clopath et al. (2010): in silico activity-dependent 
(dis)appearance of a connection occurs in terms of a competition among "strong" links in a "sea" 
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of weak synapses, which are undetectable by cehular electrophysiology. There is still no general 
agreement on whether STDP and its variants can entirely account for developmental wiring, fine- 
tuning, or remapping of microcircuits connections (j72p . However, the role of STDP for connectivity 
development has been positively discussed for systems where information is present at fast time 
scales (j73[) . as it might be relevant for the neural code of cortex and olfactory microcircuitry. In 
addition, the STDP "triplet" rule proposed by Pfister and Gcrstner (|35p . which is the central 
ingredient for this work, captures the behaviour of developmental plasticity models ([37]) widely 
used in classic studies of connectivity development (j74p . As a consequence, our results are built 
upon well-known phenomenological models, which have been shown to have excellent agreement 
with experimental data sets (|21l [35l [30)) . We have, therefore, neglected more accurate biophysical 
descriptions of SD ([ZlJ [TH] [77l [TH]) and STDP (j79l[80l), aiming at a minimalist description of neu- 
ronal excitability, synaptic transmission, and synaptic plasticity. The advantage of our approach 
is that the functional consequences of the interactions between SD and STDP could be analysed 
by standard mean-field analysis ([Ml . 

Our study concludes that the symmetry of the connections may be influenced by the internally 
generated spatio-temporal correlations of neuronal firing (|8ip . Therefore, connectivity motifs might 
not emerge exclusively from correlations governed by the external inputs pop . Both internally 
generated and externally imposed correlations are likely to play a role and therefore the results of 
Clopath et al. (2010) still hold. For instance, in a microcircuit connected by purely depressing 
synapses, external activity can still induce heterogeneity in the firing rates. If a subset of neurons 
was strongly depolarized by external selective inputs, then units would fire above the critical 
frequency, and non-random bidirectional connectivity might equally emerge (see Fig. [5]). Similarly, 
strong hyper polarization by external selective inputs, would cause neurons to fire below the critical 
frequency, even though neurons are connected by recurrent facilitating synapses (not shown).. 

Out of many possible features affecting internally generated activity, such as cellular excitability, 
architecture of long-range connections, and SD properties, here we focused on the last one. We note 
however that the dataset of Wang et al. (2006) contains additional information on cell excitability 
and its correlation to motif emergence. Cells displaying accommodating discharge patterns (i.e., ex- 
hibit spike frequency adaptation) are found to establish mostly unidirectional motifs. Instead, cells 
displaying nonaccommodating discharge patterns establish mostly reciprocal motifs (|19p . A differ- 
ential expression of spike-frequency adaptation currents results in (non)accommodating discharge 
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patterns. This is known to affect the firing rate distributions in recurrent model networks (j821 1831 
and generaUy lower the location of steady-states stable equilibria (Figs.Hp-E). For the model 
parameters employed here, we did not observe differences when repeating all our simulations with- 
out spike-frequency adaptation. Adaptation currents in the model are not strong enough to reduce 
substantially the steady-state firing rates of e.g., the facilitatory subnetworks below the critical 
firing rate for STDP, and violate our conclusions. On the contrary, as spike-frequency adaption is 
reported by Wang et at in neurons participating to unidirectional motifs, we speculate that it par- 
ticipate jointly with depressing SD in maintaining STDP in its correlational "pre-post" mode: once 
more this determines the emergence of unidirectional connectivity motifs. Non-accommodating dis- 
charge patterns would by definition not interfere with the output firing rates, in networks connected 
by facilitating synapses and bound to emerge reciprocal motifs, by the "Hebbian" mode of STDP. 

Our study addresses the emergence of motifs in the cortical pyramidal microcircuitry, but as 
preliminary data collected in the olfactory bulb become available, our framework could tentatively 
model a common principle for synaptic wiring. Long-term and short-term plasticity has been ex- 
perimentally found among olfactory mitral cells, ((5^ [2^ . and STDP was reported in the rodent 
and insect olfactory systems ([551 [5^ , and invoked at the output of mitral cells to explain decor- 
relation of sensory information (j87p . We, therefore, speculate that similar mechanisms for the 
emergence of connectivity motifs arc common in both the cortical and olfactory microcircuitry. 

Our investigation of heterogeneous network architectures revealed that STDP and SD led to 
a stable configuration of developing connections, compatible to the experimental observations. 
Under the simplified hypotheses that we adopted, the configurations of connectivity motifs and 
their relationship to the underlying neuronal collective dynamics arc stable. Our results have 
been grounded on both numerical simulation results and mean-field analysis. The latter, however, 
ignored spike-timing correlations between spike-trains. It could thereby only assess the conditions 
for the emergence of strong connections. While this is sufficient to anticipate the emergence of 
reciprocal connectivity motifs, its converse does not per sc implies the emergence of unidirectional 
motifs. In such a case, connections would for instance tend to weaken or to disappear. It is only 
through the numerical simulations of the full network of spiking neurons that we could show that 
the causal spike-timing information, if present, gives rise to unidirectional connectivity via STDP 
mechanisms. 

We also note that the major difference in SD properties, which accounts for motifs emergence, 
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is the heterogeneity of the time constant representing the short-term depression recovery Tree- In 
this respect, our results and conclusions would be qualitatively unchanged by replacing facilitating 
synaptic properties with linear, i.e., non-depressing, properties. Along these lines, we hope that 
our results might prompt new experiments on connectivity motifs. We predict that the value of 
Tree, hi a pair of connected neurons, should be inversely correlated to the occurrence frequency of 
reciprocal motifs. 

Simplifying assumptions 

Among the key simplifying hypotheses of this paper, we note that STDP was assumed to scale 
only the SD parameter G, while leaving the parameter U unaltered (j58]). This choice is consistent 
to what has been reported at distinct synapses of the central nervous system, although it might 
not be representative of all cortical areas ([5II|. Although the debate on pre- and postsynaptic 
expression of both SD and STDP is fierce, our choice of SD and STDP interaction is in part an 
arbitrary hypothesis, and it serves as a first solid ground for our conclusions. Enabling STDP 
to modify the parameter U would have partly altered in an activity-dependent manner, the SD 
profile of a synapse. This would have made isolating SD contribution more complex, and relating 
our findings to previous theoretical works (|67l l68l l69l [30l [70|) less straightforward. 

The model itself is purely phenomenological, and does not capture biophysical details, but 
rather the interaction of SD and STDP via a set of variables locally known to the synapse. It 
does however maintains the desirable compatibility with experimental data. In addition, explicit 
and systematic details on STDP at excitatory facilitatory synapses in the ferret prefrontal cortex 
are currently scarce (jl9p . While more efforts, both experimental and theoretical, should be un- 
doubtedly devoted in these directions, the hypothesis of scaling G and not U remains a simplifying 
assumption, in the view of lacking of a systematic understanding on how STDP affects all the 
parameters of the SD model (|89p . 

The presence of two classes of excitatory cortical synapses (|19p. with a distinct set of SD param- 
eters such as Tree, Tfacil, and U , prompted us to consider in our model the existence of two classes 
of excitatory neurons. The consequence of having these classes within the same microcircuitry, as 
well as their postsynaptic target preferences for connection, are crucial issues that deserve more 
experimental and theoretical efforts, but see however (IMIIUU)) . We thus explicitly neglected hetero- 
geneity in the target synaptic preference of neurons, when expressing a facilitating or depressing SD 
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profile (|66p . We instead considered homogeneous and heterogeneous nem'onal populations, where 
the facilitating or depressing SD profile was determined by the projecting neuron. Although such 
a scenario for SD seems realistic for synapses between principal neurons, it might not be accurate 
for all the synapses in neocortical microcircuits. 

The mixed microcircuitry of Fig. [7] and its preliminary analysis are an example of our attempt to 
relax these assumptions. Depending on the probability of establishing a facilitating (or depressing) 
connection, fully mixed networks were shown to behave as a continuum between the two extremes 
of the homogeneous populations studied earlier, i.e., networks of all facilitating or all depressing 
connections. These mean-field predictions were further confirmed in numerical simulations of 
large networks (not shown). In these cases, SD-STDP interactions will not necessarily lead to the 
desirable heterogeneity for the connectivity motifs, as the heterogeneity in the emerging firing rates 
does not occur. Structured local microcircuit architectures, with heterogeneous single-cell type, 
number of afferents, and firing rates are the obvious topics for a future study. We believe that 
the simple case of the heterogeneous two-subpopulation network studied in this work was a first 
necessary step towards the highlighted future directions. 

Another point that deserves some discussion is our choice for a structured, foreground periodic 
stimulation, employed in all simulations with more or less intense asynchronous background activ- 
ity. In the toy microcircuit (Fig.[T|), the cyclic stimulus was needed to demonstrate the emergence 
of a large number of asymmetric connections. It is, however, neither a strict requirement for our 
interpretation of reciprocal and unidirectional motifs formation, nor a necessary condition for their 
unbalance. This cyclic, deterministic, stimulus is an oversimplified way to invoke the temporal 
coding strategy of pop . It was chosen to unambiguously relate unidirectional pairwise motifs for- 
mation to spiking activity. Without such an input, and under the presence of strong uncorrelated 
background activity, unidirectional motifs would tend to be very sparse. We underline that our 
intention was to expose to the very same conditions networks identical in every aspects but in the 
synaptic SD, and validate the heterogeneity conditions of the connectivity matrix. 

Limits of our approach 

Our proposed mechanism for non-random pattern emergence is based on the sole interactions be- 
tween STDP and SD. Obviously, it is unlikely that these mechanisms operate independently of 
other synaptic phenomena. Homeostatic plasticity could for instance continuously rescale synap- 
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tic efficacy and make SD heterogeneities less predominant in determining connectivity motifs. In 
the lack of a priori experimental information, we chose the maximal synaptic efficacy A to share 
the exact same value for all the microcircuits we examined, in order to ensure a fair comparison. 
With all its limitations, our proposal may still provide a simple working hypothesis on one com- 
ponent underlying the emergence of connectivity, linked to short-term synaptic dynamics, along 
the same lines of the theory proposed by Clopath et al. (2010). Its validity could be challenged 
by experiments that interfere and probe the emergent firing activity, e.g., in local in vitro cultured 
microcircuits with known synaptic properties (j9ip . Our framework might be also useful for inves- 
tigating further structure-function relationships at the subcellular level, by altering the synaptic 
machinery, or by employing (future) genetically-encoded fluorescent reporters of synaptic efficacy 
and dynamics. The use of optogenetics and genetically encoded neuronal voltage- and calcium- 
sensors, may lead to experimental validation or falsification of our hypothesis, which might directly 
contribute to understand short- and long-term plasticity interactions. 

Further, the results of our simulations do not automatically capture the over-expression of 
reciprocal connectivity cortical motifs over the unidirectional ones ([7]) . This is violated in the case 
of short-term depressing synapses, where the expression of reciprocal motifs is clearly at chance 
level Fig. [SJ\, but it is however not the case for Fig.[5}3. Clopath et al. (2010) model these aspects 
by employing different neuronal coding strategies, i.e., time- or rate-codes. In our work, in order 
to isolate the heterogeneity component of SD properties, we considered a common type of spatio- 
temporal activation and ignored, for the sake of simplicity, heterogeneous neural coding strategies. 
We are however confident that including more complex network architectures, as well as more 
specific spatio-temporal correlations, as in Clopath et al. (2010), additional non-random features 
could be explained. 

In our results, we use oversimplified network architectures to illustrate that two recurrent 
excitatory networks, identical in any aspect apart of their synaptic properties, evolve different 
connectivity motifs. We simulated and studied these models analytically, and we identified the 
mechanism behind the heterogeneity as frequency-related. It would not be an obvious easy task to 
explain simultaneously all the properties of cortical connectivity, as experimentally observed (|19p . 
The choice of our homogeneous architecture, the lack of cell diversity, and the isolated evolution of 
neuronal dynamics, considerably differ from the actual cortex as well as from the variety of internal 
and external stimuli it receives during development and during lifetime. Nevertheless, the choice 
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of our protocol serves in demonstrating a potential mechanisms behind the motif formation, for 
facilitating and depressing synapses. 

We would like to emphasise that our theory refers only to one of many possible, perhaps com- 
peting, mechanisms that contribute to stereotypical motifs emergence. Alternative explanations 
and a causal demonstration of the key ideas we suggest, remain to be provided. It might be of 
interest exploring to which extent developmental changes in SD, such as the switch from depression 
into facilitation at synapses between layer 5 pyramidal neocortical neurons (|64p , occurring after 
postnatal day (P) 22, are mirrored by changes in motifs statistics. For marginal pair-wise proba- 
bility of connection. Song et al. (2005) (l7|) report no significant dependence on age, but provide no 
systematic characterization of motifs statistics beyond P20. 

It may be also possible to attempt a chronic manipulation of the firing rates of neuron (sub)populations, 
by pharmacologically altering synaptic profiles, e.g., modulating postsynaptic receptor desensitisa- 
tion, changing the presynaptic probability release, or interfering with neurotransmitter recycling. 
As future directions, more complex heterogeneous anatomical architectures and single-cell prop- 
erties should be incorporated within the same computational modeling framework. Very specific, 
non-random initial architectures, e.g., small- world and scale-free (j92p . could be explored, extend- 
ing our results towards other aspects that determine reciprocal or unidirectional motifs, possibly 
beyond the firing levels and towards, e.g., , the density of hub nodes, ranking orders, heavy tails 
in neighbours distribution, etc. 

How critical is the STDP "triplet" rule for the validity of our results? Any STDP model that is 
able to capture the high frequency effect on plasticity, as revealed by the experiments of Sjostrom 
et al. (2001), will reproduce our results if combined with short-term facilitating and depressing 
synapses. Not all STDP models are consistent with the data of Sjostrom . For instance, Froemke 
and Dan (j93p proposed an STDP model capturing pre- and postsynaptic frequency effects. This 
model predicts that in high frequencies no reversal of LTD into LTP takes place, an inversion that is 
critical for the mechanism we propose for motifs emergence. Interestingly, however, a generalization 
of the "triplet" rule that we adopt here is able to capture the experimental data of Froemke and 
Dan (2002) ([M|) . We also note that pair-based STDP models may still support the emergence 
of reciprocal connectivity motifs and replicate our results (not shown). In fact, by adjusting the 
numerical values of its parameters (e.g., = "^A^ in Eqs. S39 and S42, see Supplemental 
Information), LTP can be made prevailing over LTD above a frequency similar to the critical value 
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of the triplet-model (|72l Nevertheless, when the pair-based STDP parameters are chosen to 
match the lOHz temporal window of the triplet-STDP model employed here (compare Figs. [5]^ 
and[5J\), the frequency-dependent profile observed by Sjostrom et al. (2001) is not reproduced 
(compare Figs. [2j3 and [3lC) and SD-STDP interactions do not lead to the same results (compare 
Figs.[2]C,D and[3B;,G). This served here as a negative controlsk for our results. We also remark that 
excluding heterogeneity in spike propagation delays made our analysis simpler. It is known that 
the distribution of delays considerably affect emerging network states, besides having an obvious 
effect on STDP (|96)) . In this respect, not only STDP and SD interactions might be altered by 
propagation delays, but the collective dynamical properties of recurrent networks could also vary. 
Once more, for the sake of simplicity we chose to consider a minimalist scenario, setting the ground 
for future, more extended, investigations. 

Finally, we underline the great value of the availability of physiological information accom- 
panying anatomical connectivity. These complementary data-set contains precious statistical in- 
formation regarding the expression of microcircuit motifs, which we are starting to recognize ([TJ 
[5]). We believe that computational modeling is in this context a very powerful tool to explore 
additional hypotheses and challenge further theories. 
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A Supplementary Text 

For the sake of illustration of the standard mathematical techniques employed in the main text, 
we consider two extreme simplifications of the mean-field description of a recurrent networks with 
short-term plastic synapses (Eqs. 9-10; see also (I36p . and references therein): the predominance 
of depressing mechanisms or the predominance of facilitatory mechanisms. We analyze these two 
cases and specifically study the stability of dynamical equilibria. These simplified models and their 
analysis have been already provided in the literature, with varying degree of details. We further 
comment on the analysis of the full model of short-term synaptic plasticity and comment on the 
impact of adding recurrent inhibition to our scenario. We then consider and numerically analyse the 
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mean-field description of a mixed population, where a generic neuron may establish simultaneously 
depressing and facilitating synapses to its targets. We finally examine and partly relax the necessary 
conditions for a heterogeneous network of two interacting subpopulations to display emergence of 
connectivity motifs. For the case of a random matrix (i.e., as a null hypothesis), we provide some of 
the statistical properties of the symmetry measure we adopted to quantify the connectivity motifs 
to derive a confidence measure. In the last section, we provide model details for the pair-based, 
STDP as well as for the anti-STDP "triplet" model, discussed in the main text. 

A.l A single population with depressing synapses 

The Eqs. 9-10 of the main text can be simplified under the hypothesis of very fast recovery from 
facilitation. This is the same of assuming that Tfadi is very small and that u — U docs not vary in 
time, as a consequence. In this case, only short-term depression is modifying the synaptic efficacy, 
so that the mean-field firing rate dynamics of a large recurrent network of excitatory neurons, can 
be described by a system of two dynamical equations (|33]): 

Th=-h + AUxE + lext 

(S13) 

X ~ (I — x) /Tree ^ U X E 

These are coupled non-linear differential equations that can be analyzed by standard methods 
of dynamical systems theory (j59p . We will consider here the derivation of the equilibria for h and 
X, and indicate how their stability was assessed. By definition of equilibrium, h~0 and x = 
for those points (also called fixed points) . Substituting these conditions in Eas. lS13l we obtain two 
non-linear algebraic equations, 

h= AU X E + lext 

(S14) 

X = 1 /{\ + U Tree E) 

Using the second equation to replace x as it appears in the first, we get an implicit equation in 
the unknown h: 
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h= AU E /{l + UTrecE) + I,,t (S15) 

Eq. lSlSl can be solved numcrieally (e.g., by the Newton-Raphson method, ^7} ). given a specific 
set of parameters A, U, Tree, and lext- Ahernatively, its solution(s) can be interpreted graphicaUy 
as the intersection(s) between two functions of h, Fi{h) and F2{h), in the cartesian plane, 

Fi{h) = h 

(S16) 

F2ih) = AU E /{I + U Tree E) + h^t 

with Fi{h) is the unitary slope line (see Fig. IS1|) . Retaining the graphical interpretation, it is 
possible to appreciate intuitively the existence of the equilibrium points and their dependence 
on the parameters, as explored in Fig. ISII To this aim it is useful, prior to plotting F2{h), to 
analytically determine some of its mathematical properties, such as the asymptotic limits and 
derivatives. 

We observe that, by definition of E' = [a {h — 0)]^, F2{h) is zero for values of h lower than 9, 
h < 0. Ash ^ +00, i^2(/i) tends to an horizontal positive asymptote, occurring at (J + /ext)- 
For values of h larger than 9, the first derivative of F2{h) is always positive, indicating that the 
function is monotonically increasing. In the same range of h, the second derivative is instead always 
negative, therefore indicating that the function is convex. Moreover, the tangent line to F2{h) at 
h = 9 is the steepest of all the tangents to subsequent points (h > 6). 

F2{h) = a AU /[I + aTreeU {h - 9)f 
< (S17) 

F2{h) = - 2a^ Tree A U^ / [1 + a Tree U (h ~ 9)f 

The value of the slope of the tangent line at F2{h) aX h ~ 6 is particularly informative when 
drawing F2{h), since at any coordinates h > 9 all the other tangent lines are by definition less 
steep than it. This maximal slope {F2{9) — a A U) can then be used as a necessary condition 
for the intersections between F2{h) and the unitary slope line, at least when lext £ 

When lext < 0, there is always an intersection between F2{h) and Fi{h) a.t h — lext- This 
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is a stable equilibrium point (i.e., see below for the discussion of the stability), li a A U < 1, 
there will be for I^xt < no other intersections, since Fi{h) = h has unitary slope itself. Given 
the necessary condition a A U > 1, there exist a minimal value for those parameters, above 
which other two intersections (i.e., one stable and one unstable) with the unitary slope lines occur 
(compare Fig. ISIIA.B. and C, which were obtained for increasing values of A). 

Determining analytically such values requires imposing the condition where the unitary slope 
line becomes tangent to i^2(^)- Mathematically this can be expressed by stating that when /ea;( < 9 
there is a specific (intersection) point Hq > 9 between Fi{h) and F2{h). By definition, this point 
lays on the unitary slope line -F2(^o) = ^o, with f2(^o) = Ij i-e., F2{h) has a unitary first 
derivative at ho (see Fig. ISIB ). 

F2{ho) ^ a AU {ho-e) / {1 + a Tree U {ho - 9)) + lext = ho 

(S18) 

F2{ho) ^ a AU / [I + aTreeU {ho ~ 9)f = 1 

Simplifying the algebraic manipulations, we note the apparent similarities between the two 
equations above. We express (part of the) numerator and denominator of the right hand sides of 
each equations, by denoting the common terms as A^o and Do, respectively, as it follows 

F2{ho) = No {ho -9) /Do + hxt = ho 

(S19) 

F2{ho) = No/ Dl = I 

I. 

Thus, from the second equation we derive A^o / -Co = Do, and we substitute it in the first, 
obtaining 

ho ^ 9 + ^{9 - Ie,t) / {a U Tree) (S20) 

The above expression is of course only defined when the argument of the square root is positive, 
which is consistent with our previous hypothesis (i.e., lext < 9). One can now replace ho in the 
second equation of Eqs. IS19I and obtain the corresponding critical value of A, associated to the 
existence of such a (double) intersection (Fig. ISIB): 
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Ao = [l + y/aUTrec [0 - /ext))^ / (« U) 



(S21) 



In summary, when I^xt < there is always one (stable) equilibrium at h = I^xt and of possibly 
other two intersections (one stable and one unstable; compare Fig. ISlt \-.B. and C), depending on 
the strength of A with respect to Aq. For If-xt > 0, the situations changes and the scenario 
simplifies considerably, with only one (stable) intersection for any other choice of the parameters 
(see Fig.EljD). 

The analysis of the stability of the equilibrium points conclude our discussion. Since the 
dynamical system described Eqs. lS13l is non-linear, stability of equilibrium points must be related 
to the linearized system. The linearization is obtained for each equilibrium point by first-order 
Taylor expansion of Eqs. IS131 around that point. Let's compactly rewrite Eqs. ISlBl as 



h ^ Gi{h,x) = {-h + AUxE + hxt) /t 

X = 6*2 (ft., x) = (1 — x) /Tree + UxE 

The linearized system, around a generic equilibrium point {ho,xo), is then 



(S22) 



h w Gi{ho,xo) + dGi/dh\(^h^.xa) (h ~ ho) + dCi/ dx\(hg^xQ) {x ~ x^) 



x « G2{K,xo) + dG2/dh\(^h^^xo) (h-ho) + dG2/dx\^h„^x„} {x - xq) 



(S23) 



Assessing the stability of the above system reduces to studying the Jacobian matrix M: 



M{ho,XQ) 



( \ 

dGi/dh dGi/dx 



V 



\{ho,XQ) 



(S24) 



dG2/dh dG2/dx 
By definition, it is possible to make explicit the Jacobian matrix as 
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M{ho,xo) 



( \ 

(-1 + a A J7 xq)It {a AU {ho - 0))/t 



V 



~a U xq 



aU (ho - 0) 



(S25) 



In particular, the real part of the two eigenvalues associated to M{hQ, xq) has been analyzed for 
each equilibrium point (/io,a^o)- The eigenvalues A1.2 were computed as the roots of the following 
algebraic second order equation 



det{I - A Af(/io,xo)) = 



(S26) 



where / indicated the 2x2 identity matrix and det{) indicates the computation of the determinant 
of a square matrix. When at least one of the eigenvalue had positive real part, the equilibrium 
point was classified as unstable. When both eigenvalues had negative real parts, the equilibrium 
point was classified as stable. Nothing can be however concluded on the stability of the non-linear 
system, in the cases in which one or both eigenvalues have zero real part (and the other has negative 
real part). Stable and unstable equilibrium points have been graphically represented as circles and 
squares in Fig. [SIl respectively, as well as in Figure 4C-E. 



A. 2 A single population with non-depressing facilitating synapses 

The Eqs. 9-10 of the main text can be again simplified under the hypothesis of very fast recov- 
ery from depression r^ec- The mean-field firing rate dynamics of a single neuronal population, 
recurrently connected by short-term plastic synapses, can be rewritten as: 



T h ^ -h + A u E + le. 



[U - U)/Tfaca + U [l - u) E 



(S27) 



While this hypothesis is not as realistic as the one of the previous section, it is preparatory for 
the analysis of the full model. As for the previous case, we consider the derivation of the equilibrium 
points for h and u as well as the assessment of their stability. Substituting the definitions of 
equilibrium, h = and ii = 0, into Eqs. IS27[ we get 
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Au E + le. 



U = U {1 + E Tfacu) / {I + U E) 



(S28) 



Using the second equation to replace u as it appears in the first, one obtains an imphcit equation 



in h: 



h= AU [1 + E Tfaca) E /{l + U TfacU E) + 



(S29) 



As for Eq. IS15[ numerical methods can be used for solving Eq. IS29[ looking for values of h that 
satisfy the equivalence given a specific set of parameters A, U, Tfadi, and lext- The solution(s) 
of Eq. [829] can be also graphically interpreted as the intersection(s) in the cartesian plane of two 
functions of h, Fi{h) and i^3(/i) as defined below 



Fi{h) = h 



F3(h) ^ AU {1 + E Tfaca) E /{I + U TfacU E) + 



We observe that by definition of i? = [a (/i — 0)] , , F^ih) is zero when h < 



(S30) 



When 



h — >■ +00, the function diverges to infinity, but it can also be approximated by the straight line 
F3{h) « a A h. We also note that for values of h larger than 9, the first derivative of F3{h) 
is positive, indicating that the function is monotonically increasing. In the same range of h, the 
second derivative is also positive, therefore indicating that the function is concave. 



a A U [l + 2 a Tf^^u (h - S) + t^.^^.^ U {h - Sf 
[1 + a Tf^^ii U (h - 9)]^ 



F^{h) = 2 a2 Tfaca AU {1 - U) / [I + a Tfaca U [h ~ e)f 



(S31) 



As opposed to the previous case, the value of the slope of the tangent line at F^lh) a.t h = 6 
is not particularly relevant when drawing F^lh), since at any coordinates h > all the other 
tangent lines are by definition steeper than it. The minimal slope is F^^O) = a A U can then 
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used in combination with the asymptotic approximation F^lh) ^ a A h (i.e., the maximal slope 
is a A). It is then clear that, for < lext < ^, a sufBcient and necessary conditions for having 
always two equilibrium points (i.e., one stable at the value h = I^^t when < Ic,xt < and 
the other unstable) is represented by a A > 1 (Fig. IS2IA.B). All the considerations on how to 
assess the dynamical stability of these equilibrium points hold, and the expression of the Jacobian 
matrix Af , whose eigenvalues determine the stability, is given below: 



M{ho,uo) 



( \ 

{-l^oL Auq)It {aAU{ho - B))/t 



^ -a t/ (1 - wo) -T/l; - aU {h^ - 9) 
A. 3 Single population with short-term plastic synapses 



(S32) 



In the general case, the mean-field equations of a single neuronal population, recurrently connected 
by short-term plastic synapses, are given by 



T h = —h + A u X E + /e; 



(1 - X) /Tre 



u X E 



(S33) 



([/ - U)/Tfaca + U {I ^ u) E 



with the neuronal gain function chosen as a threshold-linear relationship between input (mean) 
current h and output firing rate E = [a [h — 9)]^ (see also Eqs. 9-10 of the main text). The 
analysis of this system, including its equilibrium points, has been already given elsewhere (|36p . 
Supporting the necessary condition on the symmetry breaking by long-term plasticities mentioned 
in the Results section of the main text, here we derive a simple observation on the analytical 
properties of these equilibrium points. According to the definition, we substitute /i = 0, i = 0, 
and u = into Eqs. IS33[ and get 
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A U X E + lext 

1/(1 + E) (S34) 

[/ (1 + ETSac^l) / (I + C/ E) 



By appropriate substitutions of x and of u into the first equation, it is possible to express it as 
h = Fi{E{h))^ an implicit equation in h: 

We observe that for h — > +oo, E{h) +oo and Fi{E) — > A Tr^J. + implying the 

existence of an horizontal asymptote. This intuitively suggests that for any choice of the other 
parameters compatible with the existence of multiple equilibrium points (i.e., intersections between 
Fji{h) and the unitary slope line), the uppermost equilibrium point (i.e., always stable) will change 
its location proportionally to A and to t^T^],, for the same choice of lext- Hence, a high value of 
Tree, as in depressing synapses, will give a lower asymptote versus a low value, as in facilitating 
synapses. 

Let's now consider two independent populations of excitatory neurons, one recurrently con- 
nected by short-term depressing synapses (i.e.. Tree > Tfaeii) and one by short-term facilitating 
synapses (i.e., TfacU > Tree) and both receiving identical non-zero external inputs lext- As for 
the previous considerations on the horizontal asymptote of F^ (E) , for an appropriate choice of A 
(i.e., large enough to have multiple equilibrium points) or for any value of lext > the firing 
rate uppermost equilibrium point of the facilitating population will always be larger than the firing 
rate uppermost equilibrium point of the depressing population. Together with the specific firing 
rate dependence of STDP, arising from the triplet-interactions, this consideration rules out that 
reciprocal motifs of short-term depressing synapses will outnumber unidirectional motifs of facil- 
itating synapses. The stability analysis for the depressing and facilitating populations (with the 
parameters used in our simulations) is provided in the main text (sec Fig. 3D-F). 

Assessing the stability of the above system reduces to linearization of the system around the 
fixed points by the use of a Taylor expansion and the study of the so called Jacobian matrix 
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AI{ho,XQ,uo) of the system, which for the sake of completeness, we report below: 



{~l + aAuoXo)/T {aAuoiho - e))/T {aAxo{ho ~ 9))/t 



-auoXo 



'T^.^l ~ auo{hQ - 9) ~axo{ho~d) 



V 



all (1 — uq) 







(S36) 



A. 4 The impact of recurrent inhibition 

We extend the description of the system given in Eqs. IS13[ to the case where recurrent inhibition 
is explicitly accounted for. The mean-field firing rate dynamics of two neuronal populations, one 
excitatory and one inhibitory, recurrently connected by short-term excitatory plastic synapses and 
by frequency-independent inhibitory synapses (as in Fig. IS3p . can be rewritten as: 



-he + U X E - A^i I + L, 



< Ti hi = -/ij + Aie U X E 



X = [1 — x) /r^ec + U X E 



(S37) 



with E 



{he - Oe)]: and / = [a^ (h. 



We consider here the derivation of the 



equilibrium points for ft-e, ^j, and x, and we indicate how their stability was assessed. Substituting 
the conditions h^ = 0, hi = 0, and i = in Eqs. IS37[ we obtain three non-linear algebraic 
equations 



he — Aee U X E — Aei I le: 



{ h, = A,eU X E (S38) 



X = I / {I + U Tree E) 



Using the second and third equations to replace x and / in the first, we get an implicit equation 
in the unknown h^: 
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We make the assumption that the synaptic couphng from the excitatory to the inhibitory pop- 
ulation is sufficiently strong Ai,, > Tree , so that short-term depression of that pathway does not 
prevent steady recruitment of inhibition at higher firing rates of the excitatory population. We also 
assume for simplicity that recurrent excitation is also sufficiently strong so that A^e > on A^i Aie . 
Under these hypotheses, when E is below a certain value, i.e., E < 9i / [U {Aie — Tree Si)], 
the above implicit equation coincides with Eq. IS15I and it can be written as if inhibition was not 
present: 



he^ 



Aee U E 
1 + U Tree E 



- Ae 



he = AeeU E /{I + U Tree E) + h^t (S40) 

Instead, above that critical value for E (i.e., and therefore for he). Eg. IS39I does not change 
formally, apart from its coefficients: 

he = AeeU E / {I + U Tree E) + /^.t (S41) 

with Aee ~ Aee ~ Oii Aei Aie and I ext = I ext + Oii 9i Agi- It is easy to verify that 

< Aee < Aee and lext > I ext • 

The critical value for the recurrent inhibitory inputs to affect the excitatory population can 
be translated into a condition on he, i.e., he > 9e + 0% / [ae U {Aie ~ Tree 0i)]- Under 
our previous hypothesis, such a critical value for he is always larger than the threshold 6e for 
the activation of the excitatory neurons. There will exist a range of activation for he above 9e, 
where the impact of inhibition is negligible. For larger activation he , inhibition kicks by step- wise 
decreasing the parameter Aee- 

All in all, the presence of recurrent inhibition in the system does not alter qualitatively the 
conclusions on the existence of equilibrium points of the mean field description. The statement 
on the separation of the uppermost equilibrium points, associated respectively to the short-term 
depressing and the short-term facilitating networks, remains true, since the horizontal asymptote 
shares the same indirect proportionality relationship with the time constant Tree of recovery from 
depression. 
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A. 5 Single population with mixed synapses 

We consider the special case of a homogeneous network of neurons, whose connections to distinct 
target postsynaptic neurons can be simultaneously short-term depressing and short-term facilitat- 
ing. For the sake of simplicity and for distinguishing this case from the mixed populations studied 
in the main text (see Fig. 4), neurons are assumed to be indistinguishable from each other. How- 
ever, every neuron has a certain probability pjj to establish a short-term depressing connection 
with its postsynaptic target. The same neuron has probability 1 — to establish a short-term 
facilitating connection to another postsynaptic neurons. Under these simplifying hypotheses, and 
by definition of conditional expected value (|98p , the mean- field firing rate dynamics of the neuronal 
population, recurrently connected by short-term plastic synapses, is 

T h -h + A [pd ud xd + (l - Pd) up xp] E + I^xt 

x'd = (1 - Xd) /TrecD " Xd E 

< UD = {Ud - UD)/TfacUo + U D 0- - Ud) E (S42) 

x'f = (1 - Xf) /TreCF " Up Xp E 

u'f = {Up - UF)/Tfacilp + Up {I - Up) E 

The cases pu = and pu = 1 have been already examined in the previous sections. For 
intermediate values oi po, an evaluation of the equilibria of Egs. IS42] has been carried out numer- 
ically, resulting in a qualitatively similar behavior to the extreme cases, with the location of the 
(stable) equilibrium points to be intermediate between those of a short-term depressing neuronal 
network and those of a short-term facilitation neuronal network. As expected, for increasing values 
of ppi the location of all equilibrium points of E{h) (if any) decreases monotonically. 
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A. 6 Emergence of motifs when initial intra- and inter-population cou- 
pling are identical 

A specific initial configuration for tlie intra- (i.e., Jff, Jdd) and intcr-population synaptic effica- 
cies (i.e., JpD, Jdf) lias been indicated in tlie main text as a necessary condition for subsequent 
"symmetry-breaking" of the spontaneously emerging firing rates, in the depressing and the fa- 
cilitating subpopulations (Fig. 4A,C). In this section, we relax such a condition and show how, 
by an appropriate external stimulation protocol, the same symmetry breaking could occur. As a 
consequence, the emerge of connectivity motifs is generally not impaired when synaptic couplings 
are set to identical values (i.e., Jpp = Jdd = Jfd = Jdf)- Other protocols and conditions might 
lead to the same configuration and here we focus on the simplest. 

We assume that the two subpopulations receive a common external input and an alternating 
pulsed stimulus component. As in a recurring traveling wave of external activity, each subpop- 
ulation is alternatively exposed to an pulsating input component, so that both facilitatory and 
depressing subpopulations are activated but never at same time. Due to intrinsic subpopulation 
properties, determined by short-term synaptic plasticities and reviewed in the previous sections, 
and as a direct consequence of the associative character of long-term plasticity (Fig. 4B,D) dis- 
cussed in the main text, this stimulation protocol leads to stronger intra-population synaptic 
coupling and weaker inter-population synaptic coupling (see Fig. ISSp . As discussed in the results 
of the main text, such a configuration is retained indefinitely, even in the absence of external 
alternating stimulation. 

In order to understand why such a stimulation protocol succeeds in developing inter- and intra- 
population coupling asymmetries, we must examine the STDP parameters (e.g., ) in its mean 
field formulation (see Eq. 11 of the man text). Let's assume that the firing rate Ef, of the 
depressing subpopulation is larger than the firing rate Ef of the facilitating subpopulation (i.e., 
say, Ed = k Ef, with k > 1). This configuration of firing rates is forced by the external input 
component, which alternates in time and across the subpopulations. Under these conditions, the 
STDP would modify intra-population synaptic coupling Jdd as follows: 

AJdd = -A^ To, El + A+ Tr, To, El (S43) 

while for the intra-population synaptic couplings, the STDP results into: 
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^Jdf 




(S44) 



AJfd 




(S45) 



It is now easy to prove that AJjjjj > AJ^p > AJpu^ hence identical initial values 
for Jdd, Jdf, and Jfd would slowly modified heterogeneously and leading to Jjjjj > J^p 
and to Jdd > Jfd- Similar considerations can be repeated when the firing rate Ep of the 
facilitating subpopulation is larger than the firing rate Ed of the depressing subpopulation (i.e., 
say, Ep = k Ed, with k > 1), concluding that all in all that the stimulation protocol would 
shape synaptic efficacies as Jd£) > Jdf, Jdd > Jfd, Jff > J_d_f, and Jf_f > Jfd, therefore 
privileging intra-population coupling to inter-population ones, as shown in Fig. IS5B . 

A. 7 Statistics of the symmetry index 

In order to quantify and describe concisely the symmetries of the emerging network connectivity 
matrix [Aij] of size N x N, we defined the following quantity (but see e.g., (|45)) . for alternative 
definitions): 



s intuitively represents the mean absolute difference between elements that are on symmetric 
positions, with respect to the diagonal of the matrix. By definition, the elements A*j are obtained 
from Aij upon first normalizing its numerical values to the maximal allowed Amax and then 
clipping them to a lower fraction z. For instance, choosing z ~ 2/3, if Aij > 2/3 Amax then 
A* J ~ Aij I Amax, and otherwise A*^ — 0. In the Eg. IS46[ M represents the number of null 
pairs \A*^, A*^ = {0, 0} as a consequence of clipping. Then, s can be rewritten in terms of an 
arithmetic average of a set of K observations of a random variable q: 



1=1 

Assuming for simplicity that each element of \Aij\ is independently drawn from a uniform 
distribution (i.e., between and Amax), the probability density distribution of q, fq{Q), can be 




N 



N 



(S46) 




(S47) 
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derived analytically (|98p . Because the arithmetic average is an unbiased estimator of the expected 
value of the random variable it samples (i.e., in this case q) ((M)) . most of the statistical properties 
of s can be immediately derived from the distribution of q. First of all, the expected value of s is 
given by the expression below 



1 - z 

1-Z2 



+ z (l + z) 



(S48) 



Deriving the expression for the variance is less straightforward and requires estimating the 
expected value of 1/K. In fact K coincides with the number of terms in the double sum of Eq. lS46l 
and it is by definition not a fixed quantity, but a realization of a binomial random variable. An 
approximated expression for the variance of s is then 



where 



1- z^ \ 6 

The validity of all the above expressions have been tested and validated numerically, directly 
estimating the average and variance of s across thousands of uniform random matrices [Aij], for 
several values of z in the range [0.1; 0.9], finding an excellent agreement. 

Finally, from the Central Limit Theorem (I98p , we can expect the density distribution of s to be 
approximately Gaussian, at least for small values of z. By the above statistical expressions, when 
studying the impact of short- and long-term synaptic plasticity in shaping microeircuit connectivity 
motifs (see Figs. 1-2,5), we expressed the significance of the observed values of s as the chance 
level, i.e., the (Gauss-distributed) probability that the observed value of s could be obtained by 
chance from a random uniform matrix. 
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A.8 Alternative STDP models 
A. 8.1 Pair-based STDP model 

By appropriately choosing A J , A2 , and setting to zero both and A'^ , earlier phenomenological 
models of pair-based STDP can be rephrased as a special case of the triplet-based model. For the 
pair-based STDP, each neuron of the network needs only two indicator variables, i.e., qi and Oi, 
instead of four. In the lack of any firing activity of the j-th neuron, those variables exponentially 
relax to zero: 

Tq, qi^ = -qij To, 61^ = -01. (S51) 

As the j-th neuron fires, the variables must be instantaneously updated. For such update rule, 
there are two distinct scenarios determining how successive pre-post or post-pre events interact 
and affect synaptic efficacy: i) all-to-all spike pairs interactions, 



gi, qi^ + I 01^ oi^ + 1 (S52) 
and ii) nearest-spike interactions. 



<7i, -> 1 01^. ^ 1. (S53) 

where the update rules do not allow accumulation of effects. Finally, when the j-th neuron spikes, 
the following updates are performed over all the indexes i: 

W^j ^ - 7] A:^ 01, (<) 

(S54) 

Wj, ^ Wj^ +77 A+ 

Following the considerations by Pfister and Gerstner pair-based models of STDP with 

all-to-all interactions must be excluded, as they do not reproduce realistic (i.e., BCM) features 
of synaptic plasticity. We then consider the triplet-based STDP model and alter some of its 
parameters as it follows: A'^ ~ A^ = 0, = 4.5 10~'^,A2 = 7.1 10^^. We also replaced 
the all-to-all spike pairs interactions by a nearest spike interaction, by modifying the update rule 
for qi and ri (and for (72 and r2, although those state variables are anyway irrelevant, upon setting 
A3 — A3 — 0). By doing so, we obtain the pair- based STDP plasticity rule matching the 
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exact same temporal window of the STDP triplet model employed here. In order to prove that 
there is correspondence in terms of the temporal window, but altered frequency-dependence, we 
subjected both the triplet-based and the pair-based models to 75 pairing events, at low frequency 
lOHz. The STDP temporal windows at such a frequency are undistinguishable from each other 
(compare Fig. 2 A with Fig. 3A). The frequency-dependence is computed across the same number 
of pairing events, imposing a pre-post or post-pre delay of 10msec. Note that for sufficiently 
large frequency of the pre-post pairs, the curves corresponding to pre-post and to post-pre become 
symmetric with respect to a horizontal line (sec Fig. 3B; i.e., around 75 — 80%, corresponding to 
long-term depression) . Such a symmetry implies that in the case of random occurrence of pre-post 
or post-pre timing, the LTP and the LTD components would cancel each other on the average. 

A.8.2 Triplet-based anti-STDP model 

This model is obtained from the triplet-based, by setting = 7.1 10^^, Aj = 6.1 10^'^, 
A2 = 0, and A2 = 3.5 10~^, by leaving unchanged the update rules for qi, q2, oi, and 02 as 
in the original triplet model, and by modifying the actual weight update equations as it follows: 
when the j-th neuron spikes, the following updates are performed for all the indexes i: 

W,j ^ W,, +T] ou{t) [A+ + A+g2,(i-e)] 

(S55) 

Wj^ Wj, - T] gi, (t) [A^ + A^ 02^ {t - e)] 

Note that this is only a tentative proposal for an anti-STDP rule, since experimental of data 
is not yet available for all induction protocols earlier employed for STDP. In particular, we ignore 
the frequency-dependency of the anti-STDP and by the new parameter set we roughly leave it 
untouched (compare Fig. 2B with Fig. 3D). 
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Figures and Figure Legends 




Figure 1. Emergence of connectivity motifs in a toy model network. Unidirectional 
(reciprocal) strong excitatory connections are indicated (A) as dashed (continuous) line 
segments, representing the topology of the network (B). Each model neuron receives periodic 
spatially alternating depolarizing current pulses, strong enough to make it fire a single action 
potential. Synapses among connected neurons display (C) short-term facilitation of postsynaptic 
potential amplitudes. Spike-Timing Dependent Plasticity (STDP) leads to strengthened 
connections and result into a largely reciprocal topology (D). Modifying the short-term plasticity 
profile into depressing (F), leads to a largely unidirectional topology shaped by STDP(G). 
Distinct motifs of strong connections arise from short- and long-term plasticities, due to distinct 
firing patterns (compare E and H), under identical external stimulation and initial connections. 
Parameters: Aij = AOOpA. 
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Figure 2. Statistics of motifs emergence in a toy model network. When decoupled from 
recurrent interactions, an isolated model synapse undergoes long-term changes depending on pre- 
and postsynaptic spike timing (A) and pairing frequency (B). Above a critical frequency (gray 
shading), spike timing no longer matters and long-term potentiation of synaptic efficacy (LTP) 
prevails on long-term depression (LTD). Panels C, D: The simulations of Fig.[T]were repeated 
2000 times, each time starting from a random initial topology. STDP progressively induced a 
persisting non-random reconfiguration of strong connections, quantified across time by a 
symmetry index (see Methods). Neurons connected only by short-term depressing synapses 
evolved strong unidirectional connections, corresponding to a low symmetry index. This is 
displayed in panel C, as an average across the 2000 simulations (left panel). Initial and final 
distributions of symmetry index values are also shown (right panel, gray and black histogram 
respectively). Neurons connected only by short-term facilitating synapses evolved instead strong 
bidirectional connections with high symmetry indexes (D). 
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Figure 3. STDP key-features for motif emergence. The pair-based STDP model, with 
temporal window shown in panel A and matching exactly Fig. [2]A., exhibits a different 
frequency-dependency (panel C) than the triplet-based STDP model (Fig. [2j3). Modifying the 
triplet-based STDP parameters to ad- hoc invert its temporal window (e.g., as in anti-STDP, 
panel B, compare to panel A), yet leaves its frequency-dependency and the LTD-reversal (gray 
shading) unchanged (panel D). Repeating the study of Fig. [2] with these two modified models, we 
find that i) the pair-based STDP fails to account for motifs emergence (panels E,G, compare to 
Fig-®, while anti-STDP succeeds (panels F,H, compare to Fig. [3]). Parameters: see the 
Supplemental Information. 



Connectivity Motifs by Synaptic Plasticity 56 




C D E 



unbalanced ext. inputs balanced ext. inputs balanced ext. inputs (zoom) 




h h h 



Figure 4. Mean-field analysis of firing rate equilibria, in homogeneous networks 
without long-term plasticity. The firing rate of Ironiogeneous recurrent networks, including 
short-term facilitating synapses (A) or depressing excitatory synapses (B), was studied by 
standard mean-field analysis. Average synaptic efficacies are indicated by Jff or by Jdd, 
correspondingly. Excitatory and inhibitory external inputs are modeled by a single term, taking 
positive, zero or negative values. A zero value corresponds to balanced excitatory/inhibitory 
inputs, while a non-zero value corresponds to unbalanced excitatory/inhibitory inputs. The 
steady-state firing rate (i.e., E, in a.u.) are the roots of the equation E{h) = h (see Egs. IIOIIT] 
and the Supplemental Information), whose graphical solution is provided (C-E), for facilitating 
(gray) or depressing (black) synapses, without long-term plasticity. Panel E is a zoomed view of 
D. (Un)stable firing rate equilibria are indicated by filled circles (squares). Networks with 
facilitating synapses fire at higher rates than networks with depressing synapses, as emphasized 
by the gray shading. 
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Figure 5. Results from numerical simulations of large recurrent networks of model 
neurons with short- and long-term plasticity. Homogeneous and heterogeneous recurrent 
networks made of 1000 Integrate-and-Fire model neurons were numerically simulated, under 
identical conditions. Panel A shows the comparison of the emergence of weak or no connectivity 
pairs (indicated as "-,-"), of unidirectional strong connectivity pairs C'^-", "D,-"), and of 
reciprocal strong connectivity pairs ("-H-", "D,D") for a homogeneous network of neurons 
connected by depressing synapses: strong unidirectional depressing connections significantly 
outnumber reciprocal depressing ones. The fractions of emerged motifs (black) is significantly 
different than the null-hypothesis (white) of random motifs occurrence. Panel B repeats this 
quantification for a homogeneous network with facilitating synapses: strong connections are only 
found on reciprocal connectivity pairs ("o", "F,F") and all emerging motifs are non-random. 
Panel C repeats the same quantification for a heterogeneous network with both short-term 
facilitating and depressing synapses. Emerging motifs display highly non-random features and 
confirm that reciprocal facilitatory motifs ("f-)'", "F,F") outnumber unidirectional facilitatory 
motifs ("—>", "F,-"), and that unidirectional depressing motifs ("—?►", "D,-") outnumber reciprocal 
depressing motifs ("o", "D,D"). Panels D-F display the steady-state firing rate distributions, 
corresponding to homogeneous depressing, homogeneous facilitating, and heterogeneous networks 
respectively. The plots confirm that heterogeneity in connectivity motifs is accompanied by 
bimodal firing rates above and below the critical frequency, represented here by a grey shading. 
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Figure 6. Mean-field simulation of a heterogeneous network with short- and 
long-term plasticity. The firing rate evolution of a heterogenous recurrent network, including 
both short-term facilitating and depressing synapses (A), was estimated by numerically solving 
the corresponding mean-field equations. The average synaptic efficacies among and across 
populations, indicated by Jff, Jdd, Jfd, and Jdf, undergo long-term modification. Panels 
B,D show the long-term changes of an isolated synapse (decoupled from recurrent interactions) 
depending on pre- and postsynaptic spike timing (i.e., tpre, tpost) and frequencies (i.e., /pre, 
fpost)- When fpre and fpost are varied independently, long-term potentiation (LTP) and 
depression (LTD) emerges as in associative Hebbian plasticity. This suggests that Jff and Jfd 
will become significantly stronger than Jfd and Jdd (C) and that such a configuration will be 
retained indefinitely. This was confirmed by simulations (E-F) plotting the temporal evolution of 
the firing rates Ef (black trace) and Ed (gray trace), and of the synaptic efficacies (F). The 
heterogeneity occurs by separation of emerging firing rates (Fig. |4]), as emphasized by the grey 
shading. Parameters: lext = 5, r = 10 msec, with initial conditions Jff = Jdd = 3, and 
Jdf ^ Jfd ^ ^ (see Supplemental Information for alternatives). 
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Figure 7. Mean-field analysis of firing rate equilibria, in homogeneous networks with 
overlapping short-term synaptic properties and no long-term plasticity. Panel A 
represents the sketch of a recurrent network where a clear segregation between subpopulations of 
depressing- or facilitating-only synapses does not occur. A neuron has a probability of 
connecting to its postsynaptic target by a depressing synapse, and 1 — pu of connecting to its 
target by a facilitating synapse. Panel B plots the location of the equilibria of the firing rate E, 
under distinct external inputs conditions and for increasing values oi po- For the same 
parameters of Fig. |4l stable equilibria move as a function of pu, taking intermediate values 
between the two extreme cases, i.e., pjj ~ and po = 1; compare to panels D-F of Figs. H) 
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Figure 8. External activity may still induce reciprocal motifs emergence in 
depressing networks. As in Figure 1, the synaptic connectivity matrix of a network of ten 
neurons was randomly initialized and pruned (A,B, i.e., pruning is indicated by the "X" 
symbols). Internally generated activity, as in Figures 1-2, contributes to the emergence of 
non-random unidirectional motifs, resulting in an asymmetric matrix W (C,D). However, if five 
units of the same network (gray circles) are externally stimulated above the STDP critical 
frequency (see the Results of the main text), a non-random connectivity emerges, featuring 
reciprocal motifs and a symmetric connectivity submatrix (E,F; upper left corner, dashed 
rectangle). The values indicated above panels A,C,E represent the symmetry index and its 
significance. 
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Symbol 


Description 


Value 


dt 


Forward Euler method integration time step 


0.1 msec 


N 


Number of simulated neurons 


10 - 1000 


Cm 


Membrane capacitance 


281 pF 


9leak 


Membrane leak conductance 


30 nS 


Eleak 


Resting membrane potential 


-70.6 mV 


-^reset 


After-spike reset potential 


-70.6 mV 


At 


Spike steepness of the exponential IF model 


2 mV 


Vg 


Spike emission threshold of the exponential IF model 


20 mV 


Vt 


Threshold voltage parameter of the exponential IF model 


-50.4 mV 




Absolute refractory period 


2 msec 


a 


Voltage dependence coefficient of the spike frequency adaptation 


4 nS 


A, 


Spike-timing dependence parameter of the spike frequency adaptation 


0.0805 nA 


Tx 


Time constant of the spike frequency adaptation 


144 msec 


'^syn 


Excitatory postsynaptic currents decay time constant 


5 msec 


Ud 


Release probability, for depressing synapses 


0.8 


Uf 


Release probability, for facilitatory synapses 


0.1 


"^rec D 


Time constant of recovery from depression, for depressing synapses 


900 msec 


'^rec F 


Time constant of recovery from depression, for facilitating synapses 


100 msec 


^facil D 


Time constant of recovery from facilitation, for depressing synapses 


100 msec 


^facil F 


Time constant of recovery from facilitation, for facilitating synapses 


900 msec 




STDP model LTD amplitude for post-pre event 


7.1 10-3 


^3 


STDP model LTD amplitude for post-pre event (triplet-term) 





At 


STDP model LTP amplitude for pre-post event 





At 


STDP model LTP amplitude for pre-post event (triplet-term) 


6.5 10-3 


Tqi 


STDP model decay time of presynaptic indicator qi 


16.8 msec 


Tq2 


STDP model decay time of presynaptic indicator (72 


101 msec 


To, 


STDP model decay time of postsynaptic indicator oi 


33.7 msec 


T02 


STDP model decay time of postsynaptic indicator 02 


114 msec 


Ai j 


Maximal synaptic efficacy 


6 ~ 12 pA 




Upper boundary for STDP dimensionless scaling factor Wij 


5 


e 


Threshold of the frequency-current response curve for mean-field models 


3 


V 


STDP plasticity rate 


1 



Table 1. Parameters employed in the simulations: STDP parameters are as in the minimal 
all-to-all triplet model described in Pfister and Gerstner (2006); short-term depression and 
facilitation parameters as in (Wang et ai, 2006); neuron parameters are as in (Clopath et ai, 
2010). 
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Figure SI. Graphical representation of the function F2{h), for different values of A 

and lext- The function F2{h) has been plotted for different values of the maximal synaptic 
efficacy A and external input lext, resulting in one intersection only (A - A ~ 3, lext ~ 0) , 
three intersections with two of them coincident with each other (B - A = 5.488, lext = 0), 
three distinct intersections (C - A = 7, lext = 0), and finally one intersection for larger 
external input (D - A ~ 3, lext ~ 4). The remaining parameters were: U = 0.8, 
Tree = 500 msec, 9 = 3, a = 1. The scripts to generate these plots and to carry out 
asymptotic analysis on the stability of the equilibrium points, see the text, are available online 
from the ModelDB (accession number 143082). 
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Figure S2. Graphical representation of the function F^^h), for different values of A. 

The function -Fa (ft.) has been plotted for different values of the maximal synaptic efficacy A and 
external input lext, resulting in one (A - A ~ 0.5, lext = 0) or two intersections (B - ^ = 3.5, 
lext ~ 0). The remaining parameters were: U — 0.1, Tfadi = 500 msec, = 3, a ~ 1. The 
scripts to generate these plots and to carry out asymptotic analysis on the stability of the 
equilibrium points, see the text, are available online at the ModclDB (accession number 143082). 




Figure S3. Mean field description for two recurrently connected populations of 
excitatory and inhibitory neurons. Excitatory neurons are recurrently connected by 
short-term plastic synapses and project to a population of inhibitory neurons. Inhibitory neurons 
project back to the excitatory population with non-plastic, linear synapses. 
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Figure S4. Graphical representation of the equivalent of functions F2{h) and F^{h), in 
the presence of inhibition. Increasing the coupling between the excitatory and the inhibitory 
population substantially bends down the curves previously analyzed as F2{h) and F4{h), lowering 
the firing rate of the equilibrium points {Aei — Aie / 10; A - I^xt = 0, B - I^xt = 5). 
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Figure S5. Formation of asymmetric intra and extra population synapses. The 

subopulations D and F of Fig. 4A receive a common input current {Ic = 12.5) and an 
alternating stimulation component {Ied, Ief), oscillating periodically between two amplitude 
levels (i.e., 2 and 0) every 1 [a.u.] of time: whenever Ied = 0, Ief = 2, and vice versa. The 
temporal evolution of the firing rate of each population is shown in A, while the corresponding 
long-lasting plastic changes of the synaptic coupling (i.e., Jffi Jdd, Jfd, Jdf) is plotted in B, 
upon initialization to the same value (i.e., 1). The same final configuration is generally obtained 
even by randomly initializing JFF^ Jdd, Jfd, and Jdf by a Gaussian distribution with mean 1 
and standard deviation 0.001, for 90 out of 100 simulation runs, demonstrating a degree of 
robustness. 
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Figure S6. External stimulation protocols and models details do not affect motifs 
emergence. The figure examines the fraction of motifs, emerging when nem'ons receive periodic 
input wave-Uke stimulation (A,B,C, exactly as in Figure 6A-C, for ease of comparison). However, 
when the periodic stimulation is omitted (D,E,F), as well as when each neuron receives instead a 
ten percent shared random background inputs with each other (G,H,I), very similar results 
emerge. The simulations of panels D,E,F, when repeated without spike-frequency adaptation 
mechanisms from each unit of the nctwork(s) (J,K,L), still give rise to the same results. 
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Figure S7. Impact of the clipping threshold parameter h on the symmetry measure 

s. The continuous trace and x markers show the expectation of the symmetry index for a 
random connectivity matrix with 20% of its elements randomly pruned, over 10000 samples of 
10 X 10 matrices, with error bars indicating the standard deviation. The simulations were 
repeated for a unitary matrix (i.e., leading to the maximum possible value for s) and for a upper 
triangular unitary matrix (i.e., leading to the minimum possible value for s). The value of 
h = 2/3 used in this work, chosen for consistence to earlier works, leads to a middle point 
between the two extremes considered and thus provide a good discriminating condition when 
using the statistics of a random matrix as a null hypothesis. 



